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PROGRAM FOR THE ANALYSIS OF TIME SERIES 


By Thomas J. Brown 

Langley Directorate, U.S. Army Air Mobility R&D Laboratory 

Christine, G. Brown and Jay C. Hardin 
Langley Research Center 

SUMMARY 

A digital computer program for the Fourier analysis of discrete time data is 
described. The program is designed to handle multiple channels of digitized data on gen- 
eral purpose computer systems. It is written, primarily, in a version of FORTRAN H 
currently in use on Control Data Corporation (CDC) 6000 series computers. Some small 
portions are written in CDC COMPASS, an assembler level code. However, functional 
descriptions of these portions are provided so that the program may be adapted for use 
on any facility possessing a FORTRAN compiler and random-access capability. 

Properly formatted digital data are windowed and analyzed by means of a fast 
Fourier transform algorithm to generate the following functions: (1) auto and/or cross 
power spectra, (2) autocorrelations and/or cross correlations, (3) Fourier coefficients, 

(4) coherence functions, (5) transfer functions, and (6) histograms. 

One of four standard data windows may be selected for application to the input data, 
and a filter, as described by the user, may be applied to the spectral data prior to output 
or generation of correlation functions. The output, as selected by the user, is written on 
a binary file for further processing or for user-designed graphics. Output is also printed 
in tabular form and/or fanfold-plot form as desired. 

The basic theory employed in the design of the program is described in sufficient 
detail to permit the user to make appropriate choices from the options available. 

INTRODUCTION 

Although there have been many computer programs written for the purpose of time- 
series analysis, each program depends upon the type of data and the objectives of the 
research. Standard Fourier series routines are useful in describing deterministic peri- 
odic functions. Through a slight generalization, a similar routine can be employed in the 
analysis of transient deterministic functions. However, when the signal is considered to 
be random, another type of analysis, based upon the concept of power spectra, must be 



utilized. Most power- spectral techniques compute the average lagged product or correla- 
tion function of the input signal, which is quite expensive in terms of time and storage. 
Further, many of these programs are essentially research tools and are inefficient for 
the analysis of large quantities of data. 

The recent advent of the fast Fourier transform algorithms has revolutionized the 
field of time-series analysis. By the proper use of these algorithms, a single program 
can be developed to handle the three types of functions discussed in the previous para- 
graph. Further, in the case of random processes, the average lagged product is no longer 
required. Thus, such a program is much more efficient than those employing the older 
technique. For this reason, digital analysis of large quantities of data becomes a 
practicality. 

This report presents a computer program for the digital analysis of random and 
deterministic time series. The program (PATS) is written in a version of FORTRAN II 
currently in use on Control Data Corporation (CDC) 6000 series machines. It employs 
the fast Fourier transform and the concept of block averaging to improve statistical vari- 
ability. It is intended for use by the practicing engineer who desires a minimum of 
involvement with the mechanics of time-series analysis. It does, however, require that 
the user possess a working knowledge of the theory of time-series analysis to obtain 
meaningful results in an optimal fashion; therefore, aspects of the theory required for 
operation of the program are discussed in this report. Only the actual equations used in 
the program are presented in the body of the report, as the basic theory is well docu- 
mented. Additional information and background may be found in references 1 to 5. How- 
ever, in the cases where the authors were unable to find satisfactory developments of the 
fundamental equations used, the necessary derivations were included as appendixes. 

Properly formatted discrete time data are analyzed by PATS through the use of a 
fast Fourier transform, from which the following functions are derived; 

(1) Auto and/or cross power spectra 

(2) Autocorrelations and/or cross correlations 

(3) Fourier coefficients 

(4) Coherence functions 

(5) Transfer functions 

(6) Histograms 

Power spectra may be filtered in the frequency domain prior to output or further 
processing with a filter of the user's description. Data may be output on a line printer in 
tabular form and/or plotted on the fanfold form, as specified by the user. In addition, all 
output is saved on binary files, which may be processed further or displayed by means of 
user-designed graphics. 
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The program was specifically designed to run economically and efficiently in a 
batch-processii^ environment from remote terminal equipment. This requirement places 
a constraint on the size of the program, which in turn limits the maximum resolution 
obtainable in the frequency domain. In some cases it was prudent to use existing sub- 
routines which are written in CDC COMPASS, an assembler level language. Functional 
descriptions of these subroutines are provided in an appendix, together with other infor- 
mation necessary to allow a user to adapt PATS to a non-CDC system with random-access 
capability. 


SYMBOLS 


®2j 

e 

E() 

f(t) 

F(cu) 

G(o;) 

> 

h(t) 

H(o)) 


Fourier coefficients 
Bernoulli number 

base value for natural logarithms, 2.7182818284 

expectation operator 

continuous function of time 

frequency of data points in amplitude bin j 

continuous function of cu, Fourier transform of f(t) 

continuous function of w, Fourier transform of f(t) defined on the 
interval T 

continuous function of w, Fourier transform of a linear system response 
continuous function of time 

continuous function of w, Fourier transform of a transfer function subject to 
input x(t) and output y(t) 


Im( ) imaginary part of a complex number 
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indices 



eq 

equivalent number of degrees of freedom 



L 

number of blocks of time data 



M 

amplitude of a square wave 



N 

number of samples of time data per block 



Nb 

number of amplitude bins 



Nt 

total number of samples of time data 



P 

period 



p(x^) 

probability density function, a function of the variable 



Pf(w) 

continuous function of w, power spectrum of the time function 

at) 


Pm 

power in the mth 1/3-octave band 



Rf(T) 

continuous function of time lag t, autocorrelation function of f(t) 


Rx(T) 

autocorrelation function (eq. (B4)) 




continuous function of time lag t, cross correlation function of 

x(t) 

and y(t) 

Re( ) 

real part of a complex number 



Sf(oj) 

continuous function of o>, power spectral density of f(t) 



Sm 

power spectral density in the mth 1/3-octave band 




continuous function of w, power spectral density of x(t) 



Sy(-) 

continuous function of w, power spectral density of y(t) 



S,y(CO) 

continuous function of co, cross power spectral density of x(t) 

and 

y(t) 
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t 


time, sec 


T time interval of length T seconds 

Ud(t),Ud(<^) data-window transform pair 

Um(t),Um(^^) Hamming data-window transform pair 

Un(t),Un(tt^) Hann data-window transform pair 

Up(t),Up(u)) Parzen data-window transform pair 

" data-window transform pair 
var( ) variance operator 

W = e-^27T/N 

Wpj data-window correction factor for correlation estimator 

Wy data-window correction factor for spectral estimator 

x(t) continuous function of time t 

X(o)) continuous function of w, Fourier transform of x(t) 

X,p(w) continuous function of w, Fourier transform of x(t) defined on the 

interval T 

y(t) continuous function of time 

Y(w) continuous function of o>, Fourier transform of y(t) 

Zj sequence of complex numbers generated from discretized time histories 

N-1 

Zj^ sequence of complex numbers related to Zj 

j=0 
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a 


(oj) 


r{) 


60 


M 


..2 


Ct) 








significance level of a distribution 
factor dependent on window chosen 

continuous function of w, coherence function of x(t) and y(t) 

incomplete gamma function 

Dirac delta function 

change 

mean value 

positive integer 

variance of the random process x 
discrete phase angle, deg 
chi-square random variable 
critical value of x^ 

expected value of x^ 

frequency, rad/sec 

discrete frequencies, rad/sec 

Nyquist frequency, rad/sec and Hz, respectively 


Mathematical notation: 


estimated quantities 


new variable 


complex conjugate operator 
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(): 

factorial operator 


ABBREVIATIONS 

DFT 

discrete Fourier transform 

FFT 

fast Fourier transform 

PATS 

program for the analysis of time series 

PSD 

power spectral density 

SFT 

slow Fourier transform 


THEORY AND EQUATIONS 
Discrete Fourier Transform 

Generalized harmonic analysis begins with the calculation of a Fourier transform, 
which assumes a definition of a transform pair. For the purposes of this program, the 
transform pair is given by 


F(cu) = 



( 1 ) 



F(o))e^“^ dcu 


( 2 ) 


When the integral in equation (1) exists, it defines a function, generally complex, known as 
the Fourier integral, or transform of f(t). The function f(t) is then known as the 
inverse Fourier transform of F(o)), and f(t) and F{o)) are said to be a transform 
pair. 

The finite Fourier transform is an approximation to equation (1) which assumes that 
f(t) is identically zero outside the region of definition. If f(t) is known on the interval 
-T/2 to T/2 continuously, then the finite Fourier transform of f(t) is given by 


t 2 ttJ_t/2 


dt 


(3) 
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When f(t) is known at N equally spaced discrete points covering the entire interval T, 
Frp(a;) may be approximated at the frequencies 


,, _ 2irk 
= 

^ N At 


(k = 0, 1, 2, . . .,N/2) 


by the discrete Fourier transform (DFT) given by 




-i2n^jk/N 


( 4 ) 


where At is the time sampling rate. This discrete Fourier transform is the basic rela- 
tion which must be evaluated in all types of digital time-series analysis. 

There are two inherent limitations in using the discrete Fourier transform as an 
estimate of the true Fourier integral. First, the finite Fourier transform assumes that 
the function for which the transform is desired is zero outside the region T. This intro- 
duces an error in resolution which is discussed in a subsequent section, "Data and 
Spectral Windows." Second, it can be shown that for N input time points, only N/2 
unique frequency points will be generated. The highest of these u>^ = 77 / At, which occurs 
when k = N/2, is called the Nyquist or folding frequency and is significant in that any 
energy present in the data with a frequency above will appear erroneously at a lower 
frequency. This phenomenon is known as aliasing and is to be avoided. Since the Nyquist 
frequency is a function of the sampling rate, aliasing may be reduced by choosing the 
sampling frequency at twice the highest frequency for which nonnegligible power occurs 
or by low pass filtering the signal at the Nyquist frequency. 

Assuming that the Nyquist frequency has been properly chosen, the discrete Fourier 
transform (eq. (4)) may be obtained from calculations of the standard relation 


N^l 

( 5 ) 

i=o 

where W = and z. is a sequence of complex numbers. To evaluate equa- 

O J 

tion (5), N^ operations are required (N multiply-add operations which must be repeated 
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N times). Implemented in this form, equations (4) and (5) are referred to herein as the 
slow Fourier transform (SFT). 

The fast Fourier transform (FFT) derives its name from its computational effi- 
ciency, which requires that N = 2*^, where N is the number of points to be transformed 
and n is an integer. For this choice of N, the number of operations is reduced from 
to 2N log 2 N, and considerable time is saved. However, the restriction that N be 
a power of 2 is often undesirable; consequently, both the FFT and the SFT are imple- 
mented in PATS and may be used interchangeably as the application requires. 

Spectral Representations 

As mentioned in the Introduction, there are several types of harmonic analysis in 
common usage. Which of these is preferred depends roughly upon whether the signal is 
steady or transient and whether it is considered random or deterministic. However, all 
these cases may be analyzed by means of the discrete Fourier transform, which was dis- 
cussed in the previous section. 

Fourier coefficients of periodic functions .- Suppose f(t) is periodic with period p. 
Then f(t) may be represented by the Fourier series 


f(t) 



cos cuj^t + Bj^ 


sin 


"k') 


( 6 ) 


where = 27rk/p are harmonics of the fundamental frequency of the signal. If N 
samples of this signal at equal intervals At are available for a total record length of 
T = N At and if T = ip, where is a positive integer, then employing these values in 
equation (5) yields 



) (k = 0, 1, 2, . . .,N/2i^) (7) 



J 


where " indicates an estimate of the required quantity. An estimate of the phase 
at frequency can also be obtained from 


®k 

4 >. = arc tan — 
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Two factors should be particularly noted in this representation: First, the total 
signal length should be a multiple of the period of the signal. If this is not the case, then 
the frequencies at which the finite transform is evaluated will not correspond to the fun- 
damental frequencies in the signal, and smearing will result. Second, since the Nyquist 
frequency occurs at the frequency o) = TxN/vp, the frequency content of the signal should 
be limited by means of a low-pass filter, because aliasing can cause significant errors in 
the estimate if the signal contains power above the Nyquist frequency. A thorough dis- 
cussion of this representation is given in appendix A. 

Amplitude spectra of transient functions . - If f(t) is a transient function, that is, it 
begins at a finite time and dies away after some time, it may be represented by the 
Fourier integral in equation (1); 


F(cu) = -i- f fWe’^^^* dt 

27t ‘^-oo 


If N samples of this signal at equal intervals At are available, a spectral estimate 
may be obtained by employing these values in equation (5). Then, by equation (4), 


^tKI = (-1)^ — Zk = 0, 1, 2, . . ., N/2) 

■‘' \ / 2ti 


( 8 ) 


where 


,, _ 27rk 
^ N At 

Power spectra of random processes . - If the signal f(t) is not considered to be a 
deterministic function, but merely one member of the ensemble which comprises a ran- 
dom process, the concept of power spectra must be employed to provide a harmonic 
representation of the function. Strictly speaking, such a representation is valid only 
when the random process may be said to be both stationary and ergodic. Briefly, this 
means that the statistics of the process are independent of time (i.e., no change in the 
mechanism of generation is present) and that each sample function is representative of 
the whole ensemble. 

When these conditions are satisfied, the Fourier integral representation given by 
equation (1) does not exist, since the function is not square integrable. However, the 
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finite Fourier transform given by equation (3) does exist and an estimate of the power 
spectral density of the random process f(t) may be obtained from 


Sf(^) - ^ 




2 


If N values of the function f(t) exist at equally spaced intervals At and these 
are employed in the standard transform (eq. (5)), then the spectral estimate becomes 



( 9 ) 


It will be shown in a later section that the factor At/47rN must be modified for other 
considerations. However, equation (9) does show the basic dependence of the spectral 
estimate on the standard transform given by equation (5). 

It should be noted that if the random process does not satisfy the condition of sta- 
tionarity, a representation in terms of power spectra is invalid and, in fact, no general 
representation of reasonable utility exists. If only the condition of ergodicity is violated, 
the power spectral representation is valid. However, many sample functions must be 
collected and an ensemble average taken over them. The reader may find a more detailed 
discussion of stationarity and ergodicity and their implications in reference 5. 

In this section, it has been indicated that the three most widely employed spectral 
representations may all be estimated from the standard transform given by equation (5). 

In the next few sections, some particular aspects of this technique will be discussed. 

Data and Spectral Windows 

One inherent limitation in techniques of spectral estimation is that the data input 
must always be finite in length. This causes a frequency smearing, or lack of resolution. 
The phenomenon may be analyzed by investigating the relation between the finite Fourier 
integral Frp(w) and the infinite Fourier integral F(aj). 

To do so, rewrite equation (3) as follows; 


F-r(w) = y ^ dt 


io)t 


( 10 ) 
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ur/2(t) 


1 1 1 t 

-T/2 T/2 

Figure 1.- The boxcar data window. 

where the "boxcar" data window, as shown in figure 1. From the Fourier 

frequency convolution theorem, equation (10) may be rewritten as 

pOO 

F((o') 

where F(w) is the true transform and the transform of Urp^gW is given by 

Thus, F,j,(ci)) is seen to be the weighted average of the values of F(w) about w = w'. 
As can be seen in figure 2, Frp(w) is an estimate of F(w), the true Fourier transform 
of f(t). Because of the duality of time-domain multiplication and frequency -domain 
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convolution, the finite transform Frp(co) at u> = is an infinite sum of contributions 
selected from F(w) by Urpy 2 (^)’ magnitude of these contributions is dependent 
upon the lobes of Urp ^2 either side of the maximum, known as side lobes. It is thus 
desirable to minimize the size of the side lobes of U,p ^2 order that F,p(u)) may 
approximate F(o)) accurately. It is worthy of note that as T — ®°, the approximation 
improves. This phenomenon may be seen by reference to equation (12), which shows that 
the main lobe narrows and the side lobes decrease with increasing T. The finite Fourier 
transform may thus be considered a smoothed approximation as seen through a window, 
in this case which is referred to as the boxcar spectral window. because of the 

characteristic square shape of its transform in the time domain. 

These data windows and spectral windows exist because of the finite length of the 
data over which the user has no control. However, there are data windows for which the 
side lobes of the corresponding spectral window are lower than for the boxcar window, 
and the averaging is thus concentrated at points nearer co\ It is therefore often advan- 
tageous to employ one of these windows. 

PATS provides three window options in addition to the boxcar. They are the Hann 
data window given by 


Un(t) = 



+ cos 


27rt 


the Hamming data window given by 


(t < -T/2)1 


(-T/2 ^ t S T/2) \ 
(t > T/2) 


0 


u 


m 


(t) = <0.54 + 0.46 cos^ 


1 " 

and the Parzen data window given by 


(t < -T/2) 
(-T/2 ^ t § T/2) > 
(t > T/2) 


1 - ^1 - 


(|t|ST/4) 


(13) 


(14) 


(15) 


V, 


(|tl>T/4) 

(|tl>T/2) 
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A thorough discussion of relative merits of the Hamming and Hann data windows 
may be found in reference 3. The Parzen window is, however, unique and warrants some 
discussion. Reference 2 will show that Parzen windows possess no negative side lobes as 
do Hamming and Hann windows. This unique feature precludes the presence of negative 
spectral estimates, which may occur when using either of the latter windows for comput- 
ing power spectra. This advantage is offset by the complexity of the window and the 
extra computing time encountered in its use. 

Because of the existence of these windows, it is necessary to use a different power 
spectral estimate from that given by equation (9) in order for the estimate to be power 
preserving. This new estimate is defined by 



(16) 


where 


Wu = ^d^(t) (17) 

is a window correction factor. The derivation of this estimate is given in appendix B. 

Variation of Estimates 

When f(t) is considered to be a random process, the power spectral estimate 
Sj(co) will be a random variable for each frequency w, because the estimate is calculated 
from a single sample function while the true spectrum is an average over the entire 
ensemble. Thus, the estimate will vary about the required value. In order to reduce this 
variation, PATS uses the concept of "block averaging." The total record length of 
points is divided into a number of blocks of length N. Then, it is shown in reference 1 
that 



where E and var indicate the ensemble expectation and variance of the random vari- 
able S£(o)), respectively, and ^ is a factor dependent upon the window chosen. 
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One simple way of assessing this variation is to assume that Sj(o;) is a chi-square 
random variable. Then, 



where kgq is the equivalent number of degrees of freedom of the random variable 
Sj(co) , and bounds on the variation may be obtained from the relation (see ref. 5) 


k SjCo;) 

< E 

^keq\2j 


r~ 1 ^ea 

[s,(o,)J 5-£2-i— 


Xk 


eq 


1-" 


(18) 


with probability 1 - a, where a is the significance level and 


\q(o;) = b 


I 


b \ eq/ eq 


The function pjxjj ] is the probability density function for a chi-square random variable 


eq 

with kgq degrees of freedom. 


The number of blocks used determines the number of degrees of freedom. For L 
sequential nonoverlapping blocks, as shown in figure 3, the equivalent number of degrees 

2Nf 

of freedom kgq is shown in reference 1 to be , that is, ^ = 1. For (2L - 1) blocks 

^ 36 Nt 

overlapping by 50 percent, as shown in figure 4, reference 1 shows that kgq that 

is, I® = ts dependent upon the window used. In PATS, kgq is calculated exactly 

through the use of equations to be found in reference 1. 


As can be seen, the number of degrees of freedom available for a fixed record length 
may be improved slightly by overlapping data blocks by 50 percent. This technique is pro- 
vided as a program option but is not recommended if the available data are of sufficient 
length to obtain the desired variance without the use of overlapping. 
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Figiire 3-- Nonoverlapping blocks. 


f (t ) 



Figure 4.- Blocks overlapping by 50 percent. 


One-Third-Octave Power Spectra and General Power Spectra 

Power spectra, as differentiated from power spectral densities (PSD), are in gen- 
eral computed by multiplying the PSD by the bandwidth of the estimate. For a sam- 
pling rate of At and block size of N points, the bandwidth of the estimate is given 
by <jo = 27t/N At. Thus, the power spectrum is computed from equation (16) as 



( 19 ) 
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One-third-octave spectra are computed by summing the contributions of the narrow-band 
spectra given by equation (16) over the 1/3-octave band. Thus, where band m is of 
width then P^, the 1/3 -octave power for band m, may be approximated by 



The 1/3 -octave power spectral density is then given by 


(20) 




2ir Aojm . 


(21) 


Estimated Cross Power Spectra 

The cross power spectral density between two signals x(t) and y(t) is estimated 
by PATS from the following equation: 



Mi 

27tWu 


z. z 


k^k 


where 


N-1 


= ^ x(j At)W^^ 


j=0 


N-1 


= J y(j At)wj^ 


j-0 


( 22 ) , 


(23) 


Block averaging is used. However, unless the coherence function, to be discussed in a 
later section, is unity, x and y are unrelated, and techniques for estimation of variance 
are not currently available. When coherence is near unity, equation (18) may be applied. 


Estimated Autocorrelation Functions 

The estimated autocorrelation function may be obtained from the PSD estimate as 
follows: 
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(24) 


Rf('T) = \ Sf(uj) do) 

Thus, the estimated autocorrelation function would be computed by applying the inverse 
FFT to the PSD estimate. This technique is appreciably faster than the method of stan- 
dard lagged products originated by Blackman and Tukey (ref. 3). 

It should be recalled, however, that since only a finite record was utilized in the 
FFT, a frequency window was introduced in the spectral estimate. As a result, the auto- 
correlation function computed from equation (24) will be distorted for large values of t. 
Thus, it is necessary to introduce a new estimate 

Rx(t) = W„ f Sjj(w)e^‘^^ do) (25) 


where 



r dt 

*^-00 


Note that the correction factor Wp^ is a function of the lag r and the data window u^(t) 
chosen. The derivation of this factor may be found in appendix C. 

An additional source of error is referred to as circular correlation error, a thor- 
ough discussion of which may be found in reference 2. Briefly, because of the periodic 
nature of the DFT, a correlation function obtained by inverting the power spectrum tacitly 
assumes that data outside the known interval are repeated periodically. This assumed 
periodicity introduces errors in the estimates for all values of lag greater than zero, as 
illustrated in figure 5. 


. N 

> 

/ ' 

/ > 

/ > 

/ 

// 

f/ 
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A 

\ 

\ 

\ 

i ■' 

kAt overlap 

_ ! 
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V. J 

Virtual data 

Real data 

-V" 

Virtual data 


Figure 5-- Illustration of circular correlation error. 
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For a lag of k At, the real and virtual data overlap at k points to introduce an 
erroneous result. A remedy for this situation is to insert zeros in the last half of the 
data block, as shown in figure 6. In the situation depicted in figure 6, the virtual data are 
set to zero; thus, a zero result is produced for the k overlapping points. PATS provides 
for zero insertion as a selectable computation option. 



/I 

r\ , 



V / 

kAt overlap 

^ i 

TaT 

/ 

N / 

Zeros 

Data 

Zeros 


Figure 6.- Illustration of zero insertion for correcting circular error. 

Spectral Filtering and Narrow-Band Correlation Functions 

Once the Nyquist frequency is chosen and a digital tape generated, a DFT analysis 
becomes somewhat inflexible. The spectra and correlation functions will contain all the 
information up to the Nyquist frequency. Often it is desirable to track a narrow band of 
frequencies, as in the case of time-space correlation studies of structures. This may be 
accomplished by digital filtering in the frequency domain. If G(u)) is the transfer func- 
tion of the desired filter and X(w) is the FFT of the input signal, the transform of the 
filtered signal Y(co) is given by Y(o)) = X(cu) G(cu). The filtered or smoothed correla- 
tion function may then be obtained by application of equations (16) and (24). 

PATS permits the user to construct a spectral filter of any general description by 
selection of a number of points on the desired response curve. A smooth function of fre- 
quency is then generated by fitting these points with a set of straight lines and quadratic 
and cubic curves in such a manner that the result is a continuous function through the first 
derivative. A thorough explanation of the technique may be found in reference 4. It 
should be noted that this option can also be employed for prewhitening and postdarkening 
should the user so desire. A thorough discussion of these filter, applications may be found 
in reference 3. 


Estimates of Cross Correlation Functions 

Cross correlation functions are computed from cross power spectral densities by 
inversion through the use of the FFT. That is. 
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(26) 



(T) = Wj^ 



Sjjy(o))e^^^ du) 


Zero insertion is optional in the program, although circular errors result if it is not 
employed. If zero insertion is employed, accurate estimates of both amplitude and phase 
are obtained with PATS. Spectral filtering for the purpose of smoothing of cross correla- 
tion fvmctions, as previously described, is also available. It should be noted that care has 
been taken so that phase is unaltered by the filtering process. 

Transfer and Coherence Functions 

The coherence function is a real-valued quantity which m"'v be estimated as 



(„). . Mr . 

Sx(o>) Sy(o)) 


(27) 


where x(t) and y(t) are the input and output of a system, as shown in figure 7. Here 
H(co) is the Fourier transform of the system response h(t). Since 


S (w) 
xy ' ' 


2 .. 

S Sjj(w) Sy(co) 


^28) 


then rxy(<^) = 1- 

2 

In the event that the system is linear, = 1 and an estimate of the transfer 
function of the linear system may then be computed by 


^(u>) 


(29) 


It should be noted that the estimate of the transfer function is valid only when the coher- 
ence is high. Also, the estimation for coherence is highly biased for small statistical 

accuracy (small number of degrees of freedom). Thus, a large number of blocks should 

2 

be averaged to get as accurate an estimate of y^^ as possible. 
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x(t) 


yit) 




Figure 7-- System response. 


Histograms 

A histogram is used to obtain estimates of the probability density function of the 
time data as follows; Consider a block of data as shown in figure 8. An array of ampli- 
tude bins is set up by PATS, as shown, and the number of points in each bin fj is counted 
from a total sample of points. The resulting histogram is then plotted. 


f(t) 



Figure 8.- Amplitude bins for histogram computation. 


The hypothesis that the process is Gaussian may also be tested by using the chi- 
square goodness -of-fit tests. To do so, the Gaussian distribution with the sample mean 
/I and sample variance is generated. The sample mean J1 is estimated from the 


Nt ^ 

original data as ^ = n and the sample variance as 

j-1 * 



The bin frequency f j is then subtracted from the expected frequency N^-Pj , where Pj 

is the Gaussian probability of the jth interval. The sum of the squares of these differ- 
ences is compared with the distribution with - sj degrees of freedom. 
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2 2 

where is the number of bins selected. The effective value of x > or is com- 
puted by 


/. .. \2 




(‘i - ^tPj) 


" jtl N,Pj 


(30) 


The probability density functions of the yf random variable are given by 


eg ^ 


(31) 


2 2 

where n is the number of degrees of freedom. The critical value of x > or Xq> lor 

the O! significance level may be found from equation (31). The normality hypothesis is 
2 2 

rejected when Xg > Xg- The decision to reject the hypothesis, however, is left to the 

2 

user, as the value of Xg is quite sensitive to Njj. A detailed description of hypothesis 
2 

testing using the x distribution may be found in reference 5. 


PROGRAM DESCRIPTION 


Operating Environment 

PATS was developed for use on the Langley Research Center CDC 6000 Operating 
System. It is written in CDC FORTRAN and uses some library subroutines written in 
COMPASS, the CDC assembler level code. The central memory requirement is 60000g 
for compiling and executing the source version presented here and SSOOOg for loading 
and executing the absolute binary version. 

Six files are used by the program: TAPEl, TAPE5=INPUT, TAPE6=OUTPUT, 
TAPE7, TAPES, and TAPE9. TAPEl is a binary file containing the input time series 
written in one of the three formats described in appendix D. TAPES is a binary-coded 
decimal (BCD) file containing the card input data in NAMELIST and FORTRAN READ 
formats. It is equivalenced to the input file. TAPE6 is a BCD file containing the output" 
to be printed. It is equivalenced to OUTPUT and is automatically printed. TAPE? is a 
binary file containing output values of all spectra, correlations, coherence, and transfer 
functions computed during execution of the program. TAPES and TAPE9 are random- 
access disk storage files used for temporary storage. TAPEl will be a magnetic-tape 
file. TAPE? may be a magnetic-tape file or a disk file copied to a tape after execution. 
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Program Specifications 


The program is written as an overlay structure with two levels. The main overlay 
sets up the COMMON storage arrays and calls the primary overlays for multiple-case 
execution. The first primary overlay reads the card input, checks for errors, prints 
informational messages, and computes the accuracy, measurement of the spectral estima- 
tors. The second primary overlay reads the input time series by blocks, computes the 
transform of each block, and stores the results on random-access disk storage. The 
third primary overlay calculates the averaged auto power spectra and autocorrelations 
and calls the plot subroutines. The fourth primary overlay calculates the average cross 
power spectra, cross correlations, transfer functions, and coherence and calls the plot 
subroutine. One large array in blank COMMON provides the temporary storage blocks 
for all overlays. Initial block addresses are assigned in the main overlay for reference 
by the primary overlay programs. Labeled COMMON blocks hold the input and control 
parameters used by all the levels. 

Appendix E contains a general flow diagram of PATS. Appendix F contains a list of 
the programs and subprograms used in PATS, with a brief description of the purpose of 
each. Appendix G contains the Langley Library subroutines used by PATS. Appendix H 
contains the source listing of PATS. 

OPERATING INSTRUCTIONS 

Deck Setups for Langley Operating System 

The following examples show deck setups for the Langley Operating System. 

Deck setup 1 .- Purpose is to fetch and execute the absolute binary version. Field 
length required is 55000g. 

JOB card 

USER card 

FETCH (A4 119, ,BINARY, ,PATS) 

NOMAP. 

LINECNT, 10000. 

1rEQUEST,TAPE 1,HY. tapeno.,ROL, 

RE WIND (T APE 1) 

SETINDF. 

PATS. 

^See Langley Computer Programing Manual for format. 
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RFL, 10000. 

DROPFIL(TAPEl) 

REWIND (TAPE 7) 

^REQUEST, TAPE 99, HY. SAVTP,RIS,your 3 initials, identification 

COPYBF(TAPE7,TAPE99) 

DROPFIL(TAPE99) 

EXIT. 

RFL,10000. 

DROPFIL(TAPEl) 

REWIND (TAPE 7) 

1rEQUEST,TAPE 99,HY. SAVTP,RIS,your 3 initials, identification 
COPYB F(TAPE 7 , TAPE 99) 

DROPFIL(TAPE99) 

27/8/9 

(Data card deck inserted here) 

36/7/8/9 

Deck setup 2 .- Purpose is to fetch, compile, and execute the source version, 
length required is 60000g. 

JOB card 

USER card 

FETCH(A4119, , SOURCE) 

RUN(S, , ,SCFILE) 

UNECNT, 10000. 

1rEQUEST,TAPE 1,HY. tape no.,ROL, 

RE WIND (TAPE 1) 

SETINDF. 

LGO. 

^See Langley Computer Programing Manual for format. 

0 

End-of-record card. 

^End-of-file card. 


Field 
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RFL, 10000. 

DROPFIL(TAPEl) 

REWIND(TAPE-7) 

1rEQUEST,TAPE 99,HY. SAVTP,RIS,your 3 initials, identification 

C0PYBF(TAPE7,TAPE99) 

DROPFIL(TAPE99) 

EXIT. 

RFL, 10000. 

DROPFIL(TAPEl) 

REWIND (TAPE 7) 

1rEQUEST,TAPE 99,HY. SAVTP,RIS,your 3 initials, identification 

C0PYBF(TAPE7,.TAPE99) 

DR0PFIL(TAPE99) 

27/8/9 

Mod card deck if any (may be a blank record) 

27/8/9 

(Data card deck inserted here) 

36/7/8/9 

Card Input Data Description 

The card input parameters are entered via FORTRAN NAMELIST and READ state- 
ments with formats for the READ statements as specified. Some parameters have 
default values noted below. Parameter types are defined as I, integer; R, real; 

A, alphanumeric. 


^See Langley Computer Programing Manual for format. 
2 

End-of-record card. 

^End-of-file card. 
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NAMELIST input format: 


FORTRAN 

name 

Default 

value 

Parameter 

type 

Description 

$INPUT 




ITFMT 

2 

I 

Code for input tape format: 

1 tape format 1 

2 tape format 2 

3 tape format 3 

(see appendix D for format descriptions) 

NFSKIP 

-0 

I 

Number of logical binary files on input 
tape to be skipped before starting exe- 
cution of this case 

NRSKIP 

0 

I 

Number of logical binary records on 
input tape to be skipped before starting 
execution of this case 

SN 


R 

Serial number of input data 

BELT AT 


R 

1/Sampling rate of input data 

STARTT 

0.0 

R 

Starting time in seconds at which pro- 
gram is to start processing data from 
input tape 

OFFSCAL 

10® 

R 

Offscale value for all channels 

NCR 


I 

Number of data channels on input tape 
(maximum value = 14) 

SCALFAC 

1. 

R 

Array of NCR values of scale factors, 
one for each channel of input data; every 
point for channel i is multiplied by 
SCALFAC (i). 

NPTOT 


I 

Total number of data points to be read 
for each channel 

NREAD 


I 

Number of data points per block to be 
read for each channel (maximum value 
is 1024 for INZERO=0, 512 for 
INZERO=l) 

LAUTOSP 

0 

I 

Array of NCR codes for computing auto 


spectra: 

1 compute auto spectrum for channel i 
0 do not compute auto spectrum for 
channel i 
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FORTRAN 

name 


Default 

value 


Parameter 

type 


Description 


lAUTOCO 0 I Array of NCH codes for computing 

autocorrelation; 

1 compute autocorrelation for channel i 

0 do not compute autocorrelation for 
channel i 

IFILTER 0 I Array of NCH codes for spectral 

filtering: 

1 filter auto and cross spectra for 
channel i 

0 do not filter spectra for channel i 
(autocorrelation and cross correlation 
will be computed from filtered spectra) 

NCROSS 0 I Number of pairs of channels to perform 

. cross functions on (maximum value = 20) 

ICROSS 0 I Array of channel numbers for cross 

functions: 

ICROSS(l,i) first channel no. for pair i 
ICROSS (2, i) second channel no. for 
pair i 

ICRSP 0 I Array of NCROSS codes for computing 

cross spectra for each pair of channels: 

1 compute cross spectrum for pair i 

0 do not compute cross spectrum for 
pair i 

ICRCOR 0 I Array of NCROSS codes for computing 

cross correlations: 

1 compute cross correlation for pair i 

0 do not compute cross correlation for 
pair i 

ITRA 0 I Array of NCROSS codes for computing 

transfer functions: 

1 compute transfer function, 

H(o;) = Sxy(w)/S^(w) 

-1 compute transfer function, 

H(w) = ^y(w)/Sy(o;) 

0 do not compute transfer function for 
pair i 
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FORTRAN 

name 

Default 

value 

Parameter 

type 

Description 

icon 

0 

I 

Array of NCROSS codes for computing 
coherence: 

1 compute coherence function for 
pair i 

0 do not compute coherence for pair i 

LAP 

0 

I 

Code for overlapping blocks of input 
data: 

1 overlap data blocks 50 percent 
0 no overlap 

IWINDOW 

1 

I 

Code for type of data window: 

0 boxcar window 

1 Hann window 

2 Hamming window 

3 Parzen window 

ITYPESP 

2 

I 

Code for type of spectral output: 

1 power spectrum 

2 power spectral density 

3 amplitude spectrum 

NPRINT 

100 

I 

Number of points to be printed from 
auto or cross spectra 

IPLOTA 

1 

I 

Code for auto spectral fanfold plots 
and/or binary tape output: 


1 no output 

2 plot and save 1/3 -octave spectra only 
(log scale) 

3 plot and save narrow-band spectra 
only (linear scale) 

4 plot and save narrow-band spectra 
only (log scale) 

5 both options 2 and 3 

6 both options 2 and 4 

•IPLOTC 0 I Code for cross spectral fanfold plots 

and/or binary tape output: 

0 no output 

1 plot and save real and imaginary 
(linear scale) against frequency 
(linear scale) 
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FORTRAN 

name 

Default 

value 

Parameter 

type 

IPLOTC 

0 

I 

FI 

0.0 

R 

F2 

20000. 

R 

LAGl 

0 

I 

LAG2 

0 

I 

PCTC 

90. 

R 

NBINS 

0 

I 

DMAX 

0. 

R 

DMIN 

0. 

R 

INZERO 

0 

I 


Description 

3 plot and save magnitude and phase 
(linear) against frequency (linear) 

5 plot magnitude (log scale) and phase 
against frequency (linear) and save 
magnitude and phase values 

Lower limit of frequency to be plotted 
on narrow-band spectra plots 

Upper limit of frequency to be plotted 
on narrow-band spectra plots 

(If FI and F2 are both zero, no narrow - 
band plots will be made) 

Lower limit of the number of time lags 
to be plotted on correlation plots 

Upper limit of the number of time lags 
to be plotted on correlation plots 

(If LAGl and LAG2 are both zero, no 
correlation plots will be made) 

Percent band to be used for calculation 
of confidence band and level of signifi- 
cance in chi-square calculation 

Number of bins to be used in 
histograms: 

0 no histograms 
(maximum value=100) 

Array of values of maximum readings 
for each channel 

Array of values of minimum readings 
for each channel 

Code for zero insertion option: 

1 insert NREAD zeros at end of input 
data block (block size is 2 x NREAD) 
(this option should be used for runs 
requesting cross correlations) 

0 no zero insertion 
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FORTRAN 

name 

Default 

value 

Parameter 

type 

Description 

NFILTP 

0 

I 

Number of points in input spectral 
filter: 

0 no spectral filter 
(maximum value=50) 

FREQF 

0.0 

R 

Array of values of frequency for spec- 
tral filter 

WGHTF 

0.0 

R 

Array of values of weights for spectral 
filter 

(user should be careful to specify points 
close together where derivative of filter 
function changes) 

$ 





Input cards after NAMELIST input: 


Card FORTRAN 

no, name 


Format 


Parameter 

type 


1 YLABEL 2A10 A 


2 TRACK 8 A 10 


A 


Description 

Array of two words (20 characters) 
to be written on each plot frame 
for case identification 

Array of NCR identification words, 
one unique word for each channel 
on input tape, eight words per 
card; more than one card may be 
needed 


Output Description 

The contents of the printed output and binary tape output of PATS are described in 
this section. 

Printed output .- The printed output consists of the following items: 

(1) Echo of input data 

(2) Informational messages about block size, type of Fourier transform to be used, 

error messages about input data 

(3) Accuracy measurement of the spectral estimators 

(4) Table of values of spectral filter calculated from input table by SPLINE 
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(5) Table of number of offscale values read for each channel 

(6) For each channel processed, 

a. Channel number, channel DD, mean and square root of variance of input 

data (MEAN and SIGMA) and the root mean square (RMS) as computed 
from the power spectral density 

b. If auto spectral output is desired, a List of NPRINT values of frequency 

and power of averaged narrow-band spectrum, and 1/3-octave band 
power spectrum and power spectral density 

c. If autocorrelation is desired, a list of time and Rx, with time lags from 0 

to N/2 (N=block size) 

(7) If histograms are requested, a fanfold plot of bin number against counts, a list of 

values of occurrences, and a goodness -of- fit test calculation 

(8) For each pair of channels processed, 

a. If cross spectral output is desired, a list of NPRINT values of frequency, 

real and irhaginary parts, and magnitude and phase of complex narrow- 
band power spectrum 

b. If cross correlation is desired, a list of time lag and Rxy with time- 

lag values from -N/2 to N/2 (N=block size) 

c. If coherence is desired, a list of NPRINT values of frequency and 

coherence 

d. If transfer function is desired, a list of NPRINT values of frequency and 

transfer function 

(9) Fanfold output - plots of every function computed for which plots are specified 

will appear in the printed output immediately following the listed values; the 
plots are limited to 256 points each to conserve line count; the first 256 points 
between FI and F2 of each spectrum are plotted; points between LAGl and 
LAG2 of correlations are plotted, skipping intermediate points to reduce the 
number of plotted points to less than 256 

Binary tape output .- Every function computed is written onto file TAPE? when it is 
computed. All calculated values are written. One record is created on the file in the 
following format; 

Word Type 

1 to 6 A 

7 I 


Description 

Label describing function and channels 
Number of points in output function, NP 
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Word 

Type 

Description 

8 

R 

First value of independent variable 

9 

R 

First value of dependent variable 

10 

R 

Second value of independent variable 

11 

R 

Second value of dependent variable 


NP + 7 R NPth value of independent variable 

NP + 8 R NPth value of dependent variable 

In printed output, a message is written noting the record number and descriptive label for 
the function. 


Restrictions and Limitations 

The restrictions and limitations for use of PATS are as follows: 

(1) The binary input tape must be positioned at the beginning of the data to be pro- 
cessed before the program starts reading data for the case. This may be accomplished 
by using control cards before execution and by assigning the correct nonzero values to 
NFSKIP and NRSKIP for each input case. For tape format 1, it should be positioned at a 
record with the desired serial number in the second word (may be after the first record). 
For tape format 2, it should be positioned at the ID record with the desired serial number 
in the eighth word. For tape format 3, it should be at the record of the file containing the 
desired serial number. If this condition is not met, a message will be printed and execu- 
tion stopped. When both NFSKIP and NRSKIP have nonzero values, NFSKIP files are 
skipped first. For tape format 2 the ID record is checked, the next two records are read, 
then NRSKIP records are skipped. For tape format 1, no records are read before skipping 
NRSKIP records. 

(2) The version of the program being presented has a maximum block size of 1024. 
The program storage requirements are SSOOOg for the absolute binary version and 60000g 
for the source version. To change this Umit, NMAX must be assigned the desired value 
in the program MAIN and the dimension of CMAIN changed accordingly. 

(3) The number of block averages and the number of individual channels to be pro- 
cessed are restricted by NBCMAX. The product must be less than or equal to 800. To 
change this limit, assign the desired value to NBCMAX in MAIN and change the dimension 
of KNDEX accordingly. 
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(4) The number of data channels on the input tape is limited to 14. To change this 
Umit, change the dimensions of ail variables dimensioned 14 in COMMON blocks BLK2, 
BLK5, and BLK8 in all overlays; assign the correct value to NCHMAX in MAIN; and make 
CMAIN dimension the larger of 4NMAX+6 or 2NMAX+6+64NCHMAX+512. 

(5) When the amplitude spectrum option is selected, no other functions will be 
calculated. 

(6) Filter input function is restricted to 50 points. To change this limit, change the 
dimensions of FREQF and WGHTF in COMMON block BLK9 in all overlays and change 
value in test in READIN (two statements after statement number 113). 

Error Messages and Remedies 
The error messages and suggested remedies are as follows: 

(1) NCR GREATER THAN NCHMAX, PROGRAM WILL NOT READ TAPE COR- 
RECTLY. JOB TERMINATED. 

To correct, see item 4 in "Restrictions and Limitations." 

(2) YOU MAY HAVE CIRCULAR ERROR IN YOUR CORRELATIONS BECAUSE YOU 
HAVE NOT ASKED FOR ZERO INSERTION. 

Job will continue. To correct, change INZERO to 1 and rerun. 

(3) BLOCK SIZE TOO LARGE FOR DIMENSIONS PROVIDED. 

Job terminated. To correct, see item 2 in "Restrictions and Limitations." 

(4) NO 50 PERCENT OVERLAP ON ZERO INSERTION RUNS. 

Input value of LAP will be altered to zero and job will continue. 

(5) INPUT INDICATES NO CHANNELS TO BE PROCESSED. 

Job terminated. Check input data and rerun. PATS resets all computing options other 
than AUTOSP to zero when ITYPESP=3. See item 5 in "Restrictions and Limitations." 

(6) NCHP*NBLK GREATER THAN NBCMAX. 

Execution ended. To correct, see item 3 in "Restrictions and Limitations." 

(7) INPUT ERROR, NFILTP GT 50. 

Execution ended. To correct, see item 6 in "Restrictions and Limitations." 

(8) NCH GT 100 NOT ALLOWED. 

Execution ended. No correction of program possible. 

(9) TAPE NOT POSITIONED AT ID RECORD FOR DESIRED SN. 

Execution ended. Correct input deck to position tape correctly and rerun. 

(10) TAPE NOT POSITIONED AT DESIRED SN. 

Execution ended. Correct input deck to position tape correctly and rerun. 
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(11) STARTING TIME NOT FOUND ON FILE. 

Execution ended. Check tape for correct times. 

(12) END OF FILE ENCOUNTERED BEFORE NPTOT POINTS READ. 

Check tape for number of points after STARTT; correct input data and rerun. 

CONCLUDING REMARKS 

This paper has presented a general purpose digital computer program for the har- 
monic analysis of multiple channels of time-history data. The program is written pri- 
marily in CDC FORTRAN and employs the technique of the fast Fourier transform. A 
complete program listing with descriptions of necessary subroutines is included so that 
the program may be adapted to any facility. In addition, the philosophy and theory 
employed by the program are discussed so that the user may make appropriate choices 
among the options available. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., April 8, 1974. 
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APPENDIX A 

FINITE FOURIER TRANSFORM OF A PERIODIC SIGNAL 

Suppose that f(t) is a signal with period p. Then, it can be represented by the 
Fourier series 


0 V 

f(t) = -^ + / (An cos + Bn sin o^n*) 


(Al) 


n=l 


O O >fr 

where ajjj = — ^ are harmonics of the fundamental radian frequency ^ of 

nal. Further, suppose that N samples of this signal at equal intervals At are avail- 
able for a total record length of T = N At. 

The finite Fourier transform of this signal is given by 


N-1 


Ik = 2 f(j At)e-^27Tjk/N 

j=0 


(k = 0, 1, 2, . . ., N/2) 


0 

at the frequencies 0)^ = These frequencies will correspond to the harmonic fre- 

quencies of the Fourier series if and only if T = tp, where is a positive integer. In 
this case, the mth harmonic will be equal to the kth frequency at which the finite transform 
is evaluated when k = t^m. 

Now, when t = tp, it can be shown that the finite Fourier transform of equation (Al) 
is given by 


Xk = — N 6(k) + N ^ — — — 6(k - im) + N ^ ^ [^6(k - wi + ^N) + 6(k + im - ZN)] 

n=l n=l i=l 


OO OO 


n=l ^ 1=1 


OO OO 


n=l 1=1 


6(k - im + IN) - 5(k + m - IN) 


where 6(j) is the Dirac delta function. 
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APPENDIX A 


The most interestii^ of these transforms are those which correspond to harmonics 
of the fundamental period, that is, k = vm. For this case. 



The summation terms in equation (A2) involve aliasing of power from higher frequencies. 
Note that the aliasing depends upon N/f, the number of points per fundamental period of 
the signal. Since the Nyquist frequency occurs when v = N/2, the aliasing may be 
removed by filtering the signal above fj, = 1/2 At = N/2ip. Assume that aliasing has 
been removed, then 


and 


"wn 


iJ 


*^m 


+ B 


m 


N 


- 

Am 


(A3) 


As an example, consider the square wave of amplitude M and period p. If the 
Fourier series representation of this signal is considered, it can be shown that 


and 



cos Wj^t dt = 0 




sin ojjjt dt = 


\4M/nn 


(n odd) 
(n even) 


Thus, the square wave admits the Fourier series representation 


OO 



n-1 


1 

2n - 1 


sin 


27T(2n - l)t 
P 


(A4) 
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APPENDIX A 


The magnitude and phase of this representation are given by 


and 


A ^ + B ^ 
^n + ^n 


r 

0 

(n even) 

[2M 

(n odd) 

niT 


. = tan'^ = tan'^(°°) = ^ 


To obtain the finite Fourier series representation for this function, assume, without 
loss of generality, that = 1 and N is an even number. Then, since equation (A4) 
yields f(0) = f(p/2) = 0, the finite Fourier transform becomes 


N-1 


Xk = ^ f(j 
J=0 


N-1 


= M ^ (e-i2-VN)' . Y (e-i^’Tk/N)' 

J=N.l 

2 


1 - e 


•ivk 


[I . e-i27Tk/N 


- 1 - M 


1 - 1 - e 




. -i27rk/N , -i27rk/N 
1 - e 1 - e 


= M 


(l - _ ^-ivk) 

-i27Tk/N 


1 - e 




-2Mi cot 


7Tk 

N 


(k even) 
(k odd) 
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Thus, 


X,, 


N 


-i2M £k 


N 


N 


(k even) 
(k odd) 


From equations (A2) and (A3), it can be seen that when k is even. 



where B 2 j is a Bernoulli number. The summation which appears in this equation is the 
aliased term which arises because the signal was not low pass filtered. 
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APPENDIX B 


SPECTRAL ESTIMATION THROUGH USE OF THE FINITE FOURIER TRANSFORM 
Let x(t) be an arbitrary stationary random .process and define 

where Uj^(t) is any data window which is zero for |tj>T/2. Then, the power spectral 
estimate becomes 


SxM=f|xr(‘^) 


' f ^2 I. “cl(*2) *(‘l) K‘2) 


-lOJ 


(U-t2) 


dtr 


(B2) 


Taking the ensemble expectation of this quantity yields 


E 


Sx<") 





(B3) 


where 


R^(t) = f Sx(w')e^^''^dw' 


(B4) 


is the autocorrelation of the random process x(t). Employing this relation in equa- 
tion (B3) gives 


p°o .. 1 r°° , \ 1 r°° -i(aj'-a))to 
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Now, define 




(B6) 


where U^(a>) will be real and even if u^(t) is even. Then equation (B5) becomes 


e[4m] = f J ^ " T I CO 

Thus, the spectral estimate obtained in this way is a smoothed approximation to the actual 
spectrum as seen through the spectral window characterized by squaring the Fourier 
transform of the data window. 

In order for this estimate to be power preserving, it is necessary for the integral 
of the mean estimate to be equal to the total power. Integrating equation (B3) yields 


£ ? 17 £ "‘l £ “d(‘l) “aW ’*x(tl - *2) at2 J_" do, (B8) 

and since 


X 1 


do) 


(B9) 


equation (B8) becomes 


£E[w]d.=i^( 

^00 pOO 

l) - tg^ 

6(t, - 1 

^ 2 ) 

_ 7T 1 1 
T 27Tv- 

..CO 

1 w 
’-00 

Rx(0) dti 




Rx(0) 

2/ 

l)% 




2T 

L “a (‘ 





(BIO) 
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Since the total power in the signal is given by Rx(0), define a new estimate 
by the following equation: 



E ^(cij) dct) = Rx(0) 


(Bll) 


Clearly, this estimate is related to the old estimate S (o)) by 


Wu 


(B12) 


where 




(B13) 


is the window correction factor. Thus, since the estimate is given by equation (9) 

as 




At |2 
47tnI 


the desired spectral estimate Sjj(oj) is obtained: 




= 1M1| 

27rWui 


(B14) 


The window correction factors for the various data windows are as follows: 
For the boxcar window, 

Wu = T 


(B15) 


For the Hann window, 
W 

Wu - 8 


(B16) 
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For the Hamming window, 


W„ = T(0.3974 + 


0.993 6 
rr 


And for the Parzen window, 


W„ =T 


151 

560 


(B17) 


(B18) 
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AUTOCORRELATION ESTIMATION FROM ESTIMATED POWER SPECTRAL DENSITY 

The autocorrelation of a function of time x(t) would normally be estimated as the 
inverse Fourier transform of the estimated power spectrum; that is, 


R, 


-(t) = r Sjj(o;)e^“'^ do) 


Thus, 


E 


^(t) = f E Sx(w) du) 
m J J . 00 L 


(Cl) 


(C2) 


Now, it can be shown from the equations of appendix B that 


e[§x(w)] = :— J ^ Sx(w') dco' 

Thus, equation (C2) becomes 


E 


Rx(^) 


= ^ r dw'e^^^”' Sx(co') f dw Uh^(u)’ - 


Now, by setting Wq = w' - w, equation (C3) yields 


E 


Rx(^) 


= ^ r, r. '*“0 ^d'(“o)' 




(C3) 




-1^0^ 


(C4) 


Recall the definition (from eq. (B6)) 
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Further, since u^(t) is real and even for all windows considered in this report, U^(w) 
is real also. Thus, 


= dtf dfu,(t)u,(t.)a-‘“(‘-*’) 

4^2 <^-00 ^-CO 


and 


r 9 1 n P 1 P 


-ict>Q(t-t'+r) 


(C5) 


However, 


277 •J-oo 


Thus, equation (C5) becomes 




1 r°° 

'2j]_„‘““d(‘)“d<*+") 


and 


E 


\ Ucj(t) Ud(t + T) dt 

R^(t) I = — ” Rx(t) 


L“> 


dt 


(C6) 


Therefore, in order to have an unbiased estimate of the autocorrelation, it is necessary 
to define the new estimate (eq. (25)): 


^(T) = W„ r do; 
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Where. 


Wr = 


r" u/(t) dt 

*^-oo ° 


poo 


is again a window correction factor. 
Note that for the boxcar window. 


and 


UT/2W - 

C "t/ 2^W dt = T 


((t I ST/2) 
(otherwise) 


r 


"T/2“> “t/2<‘ * = 




Therefore, 



(|t|<T) 


(|t|<T) 

(otherwise) 




0 


(otherwise) 
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BINARY INPUT TAPE FORMATS 


Tape Format 1 

Data digitized by analog-to-digital conversion equipment at Langley is edited, 
reduced to engineering units, and put on a digital computer tape. The computer program 
which does this is called the Adtran Quantity Pass and its standard output data tape is 
called Adtran Output Tape. The CDC Adtran Output Tape is explicitly blocked and the 
actual end-of-file mark is used to indicate the end of writing on tape. 

All tapes will be written in the binary parity using the standard CDC FORTRAN 2.0 
input/output statements. There will be between 11 and 110 FORTRAN 2.0 logical words 
per frame. These frames will be blocked into larger physical records. A file of data 
will be completely defined by serial number. New serial numbers will always begin in a 
new physical record. If a physical record is not complete, it will be filled with 999999 
(six 9's). The end of writing on the tape will be indicated by an end-of-file mark. The 
frame format is as follows: 


Word 

1 


2 

3 

4 


5 

6 

7 

8 


9 


Type 

Floating 


Contents and Description 

The number of channels of data in this 
frame; less than 40 for continuous data 
and less than 100 for commutated data 

Serial number; the input card format for 
serial number should be 6 digits wide 

Words 3 and 4 are the primary engineer- 
ing identification, for example, test and 
run; they would be represented on input 
card formats by no more than 6 digits 
apiece 

Words 5 and 6 are additional engineering 
identification 

Words 7 and 8 are Greenwich Mean Time 
and are used only for telemetry data; for 
ground facilities, word 8 may be ground 
facilities 
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Word 


Type 

Floating 


Contents and Description 

Elapsed time in seconds; processing will 
be controlled by elapsed time within a 
file; the increments in elapsed time may 
not be constant 

Data channel 1 

Data channel 2 

Data channel 3 


N + 10 


Data channel N, where N is the number 
given in logical word 1 


The relationship between frames and records is shown below. 


Number of channels, 
N 


1 ^ N £ 10 
10 < N ^ 20 
20 < N ^ 30 
30 < N ^ 40 
*30 < N S 100 


Commutated, 


Words per 
frame 


Frames per 
record 



As an example, to read a 12 -channel frame, a physical record of 510 words is read. 
The time of the first frame is in word 10, the time of the second frame is in word 40, . . ., 
the time of the 17th frame is in word 490. 


Tape Format 2 

The tape is a FORTRAN written, binary -parity, multifile tape with a flexible yet 
efficient format. Each file contains four basic record types (ID, NAMES, UNITS, and 
DATA) and consists of a continuous unique test (or run). The ID record contains non- 
repetitive information such as run or test number, date, time bias, and record blocking 
factors. The NAMES and UNITS records contain data channel names and engineering 
units, respectively. The DATA records themselves contain the engineering data. In 
addition, each record begins with a KEY word denoting the record type followed by a word 
containing the record size. Thus, all information necessary to operate on any file is 


47 










APPENDIX D 


available within the first four records of the file. The formats for the records in each 


file are as follows: 
Record 1 ID Record 




Word 

Name 

Format 


Description 

1 

KEY 

A 

ID 



2 

NN 

I 

Number of remaining words in the record = 19 

3 

IWD 

I 

Number of words of unblocked data in a data 
record 

4 

KCH 

I 

Number of words of blocked data in a data 
record 

5 

NFR 

I 

Number of frames in a data record (blocking 
factor) 

6 

m(i) 

A 

Name for first ID parameter = SERIAL 

7 

m(2) 

A 

UNITS for first ID parameter = NUMBER 

8 

ID(3) 

F 

First ID paranieter = the serial number 

9 

ID(4) 

A 

NAME(2) second ID parameter = TEST 

10 

ID(5) 

A 

UNITS(2) second ED parameter = NUMBER 

11 

ID(6) 

F 

Second parameter = the test number 

12 

ID(7) 

A 

NAME (3) = DATE 

13 

ID(8) 

A 

UNITS(3) = DAYS, YR-MONTH-DAY or 
UNKNOWN 

14 

ID(9) 

F 

PARAMETER(3) = YEAR X 10000 + MONTH 
X 100 + DAY 

15 

ID (10) 

A 

NAME (4) = BIAS 

16 

ID(ll) 

A 

UNITS(4) = SECONDS 

17 

ID(12) 

F 

PARAMETER(4) = GMT time bias 

18 

m(13) 

A 

NAME (5) = ENGR ID 

19 

ID(14) 

*I 

UNITS(5) = 2 

20-21 

ID(15-16) 

A 

PARAMETER(S) = Engineering identification 
(two words) 


*When the UNITS word for a parameter contains an integer less than 12, the param- 
eter is defined to be alphanumeric data of that many words in length. 
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Record 2 NAMES Record 



Word 

Name 

Format 

Description 

1 

KEY 

A 

NAMES 

2 

NN 

I 

Number of unblocked parameters plus num 
ber of blocked parameters 

3 NN + 2 

NAMES 

A 

Names for unblocked data parameters fol- 
lowed by names for blocked data 


Each parameter including time will have a name. 


Record 3 UNITS Record 



Word 

Name 

Format 

Description 

1 

KEY 

A 

UNITS 

2 

NN 

I 

Number of unblocked parameters plus num- 
ber of blocked parameters 

3 NN + 2 

UNITS 

A 

Units for unblocked data followed by UNITS 
for blocked data 


The UNITS are not always necessary and will sometimes be blank. 


Record 4 through EOF, DATA Records 


- 

Word 

Name 

Format 

Description 

1 

KEY 

A 

Data 

2 

NN 

I 

Number of remaining words on 
the record (IWD + KCH * NFR) 

3 

XDATA(l) 

F 

First word of IWD words of 
UNBLOCKED data (FRAME 
COUNT, e.g.) 

4 + IWD ZDATA 

ZDATA(I,J) 

F 

Blocked data I parameter, 
J frames 


The data records are optimally packed to approach, but not exceed, 512 words per record. 
The record size is determined as follows; 

SIZE = NFR * KCH + IWD + 2 

where 

NFR = (510 - IWD)/KCH 
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NFR is the blocking factor (integer) 

IWD is the number of nonrepeated words in the record 
KCH is the number of data channels 

For example, a test with 9 recorded channels and only one word of unblocked data per 
record would have 507 words in each data record as follows: 

KCH = 9 
IWD = 1 

NFR = (510 - l)/9 = 56 
Therefore, 

SIZE = 56 * 9 + 1 + 2 = 507 


Tape Format 3 

The binary tape is written by using subroutine RECOUT. The data passed to 
RECOUT at each time point are 


Word 

1 

2 

3 

4 


Contents 
Serial number 
Time 

Data channel 1 
Data channel 2 


NCH + 2 Data channel NCH 

All words are in the floating-point mode. 
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PROGRAMS AND SUBPROGRAMS USED BY PATS 


The programs and subprograms written specifically for PATS are given in the fol- 
lowing list with a brief description of the purpose of each. 


MAIN 

PLOTNB 

FANFOLD 

FOURT 

READIN 

CSQ 

FUNC 

BLOCKS 

READTPE 


Sets storage array dimensions, sets up random-access files 8 and 9, 
calls overlays for input, computations, and output 

Sets up arrays for plotting narrow-band spectra on fanfold 
References FANFOLD 

Plots an array in printed output with heading, max, min, and scale; the 
ordinate is across the page; the abscissa (index number in the array) is 
down the page, one point per line; up to 256 points per plot may be plotted 

Computes the Cooley-Tukey fast Fourier transform for an array of com- 
plex numbers; the number of points is arbitrary, although the subroutine 
operates much faster on powers of 2 

Reads NAMELIST and FORTRAN READ input, checks for input errors, 
and prints informational messages, including accuracy measurement of 
spectral estimators 
References CSQ 

Computes value of chi-square for given level of significance and number 
of degrees of freedom 
References ITR2, FUNC 

Function subprogram used by CSQ to evaluate the chi-square probability 
function 

Calls subroutines to read data from binary input tape and perform 
Fourier transforms for given number of blocks of data 
References READTPE, TRAN 

Reads one block of data from binary input tape for selected channels and 
stores the data on random-access file 9; if overlap option is selected, 
each block after first takes the last half of the previous block and fills 
the rest of the block with new data 
References WRITMS, READMS, RECIN 
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TRAN 


HANNING 

HAMMING 

PARZEN 

AUTOSP 


AUTO 


NORMAL 

PFUN 

PLOTB 

BANDS 

BNDSUM 

CROSSP 


For each channel of data selected, this subroutine reads one block of 
input data from random- access file 9, counts occurrences for histograms, 
windows the data, extends the block with zeros if zero insertion option is 
selected, performs Fourier transform, and stores the results on random- 
access file 8 

References READMS, WRITMS, HANNING, HAMMING, PARZEN, FOURT 
Weights the input array by the Hann data window 
Weights the input array by the Hamming data window 
Weights the input array by the Parzen data window 

Sets up storage arrays for subroutine AUTO and calls for histograms if 
selected; calls SPLINE to evaluate the spectral filter weighting function 
References AUTO, NORMAL, SPLINE 

For each selected channel of input data, computes the mean and variance 
of analyzed data (including overlap if used), reads all transforms for this 
channel from random-access file 8 and averages the amplitude spectra or 
PSD, applies the spectral filter, prints results, and calls selected plot 
routine; auto PSD is stored on random-access file 9; 1/3-octave spectra 
are calculated from the narrow-band spectra and printed; if autocorrela- 
tion is selected, the inverse transform of the auto PSD is performed and 
the result printed and plotted on fanfold 

References READMS, WRITMS, PLOTNB, BANDS, PLOTB, FOURT, 
ASCALE, FANFOLD 

For each selected channel of input data, calls FANFOLD to plot histogram 
data, calculates chi-square for goodness-of-fit test, and prints the results 
References FANFOLD, PFUN 

Function used by NORMAL to calculate probability density function of a 
normally distributed random variable 

Calls FANFOLD to plot 1/3 -octave spectrum 
References FANFOLD 

Integrates narrow-band spectrum for 1/3-octave power spectrum 
References BNDSUM 

Computes sum of given array of complex numbers 

Sets up storage arrays for CROSS 

References CROSS 
✓ 
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CROSS 


SPUNE 


For each pair of channels: reads transforms for both channels for all 
blocks from random-access file 8, averages the products for cross PSD, 
prints the desired results, and calls for selected plots; the cross correla- 
tion is computed from the inverse transform of the cross PSD, and the 
results are printed and plotted on fanfold; coherence and transfer function 
are calculated from the auto PSD's stored on random-access file 9 and 
the results are printed and plotted on fanfold 

References READMS, PLOTNB, ENCODE, FOURT, ASCALE, FANFOLD 

Fits a smooth curve to a set of input data points and evaluates the func- 
tion at evenly incremented intervals over a given range 
References SIMEQ 
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LANGLEY UBRARY SUBROUTINES 


The Langley Library subroutines used by PATS are ASCALE, GAMMF, ITR2, 
OPENMS, READMS, RECIN, SIMEQ, and WRITMS. The subroutine RECOUT is not used 
by PATS but must be used separately to generate input data in tape format 3. Usage 
descriptions of all these subroutines are given in this appendix. 

Subroutine ASCALE , 

Language : FORTRAN 

Purpose : To compute a scaling factor for an array of numbers to be plotted over a certain area and find 
the minimum data value within the array. 

Use : CALL ASCALE(ARRAY,S,N,K,DV), where 

ARRAY Name of the array containing the floating-point values to be scaled 

S Length (floating-point inches) over which the data are to be plotted (usually the length of 

one of the axes) 

Number of data values in ARRAY from which points are to be plotted in accordance with K 

Interleave factor which specifies the sequence in which data are stored: 

1 indicates that values are stored sequentially 

2 Indicates that values are stored in every other location in the array 

Number of divisions per inch of the plotting paper to be used (should be 10.0, 20.0, 25.0, 
or 25.4) ^ 

I 

Restrictions : The array must be dimensioned to include storage space for two extra elements per inter- 
leave factor. For example: N = 100, K = 1, DIMENSION ARRAY (102); N = 75, K = 3, DIMENSION 
ARRAY (231). 

Method : This routine scans the elements in the array to find the minimum and maximum. ASCALE com- 
putes an adjusted minimum (origin value) and stores it in ARRAY((N*K)+1) and computes a scale factor 
and stores it in ARRAY((N*K)+1+K). The scale factor will be a power of 10 x (2,4,5, or 10). The data 
in the array may be scaled to floating-point inches by using a formula similar to the following: 

SV = (AE-MV)/SF, where SV is the scaled value, AE is the present value of array element, MV is 
either the minimum value or the value desired at the origin, and SF is the scale factor computed by 
the subroutine. 


N 

K 

DV 


Storage : 262g locations for the CDC 6000 series. 

Subprograms used : ALOG, ALOGIO. 

Other coding information : Example: DIMENSION ORD(102),ABS(204);CALL ASCALE(ORD,10.,100,1,10.); 
CALL ASCALE(ABS, 15., 100,2, 10.). 


Subroutine date : September 3, 1970. 
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Function GAMMF 


Language : FORTRAN 

Purpose : To compute the incomplete gamma function 
F(A,X) = C e-f^n^-^duL 

If X = 0, then the complete gamma function is obtained. 

Use : Y = GAMMF(A,X), where GAMMF(A,X) is defined as the integral from X to «> of exp(-p) times 4 
to the (A-l)th power dp. , _ 

Restrictions : X S 0; when X = 0, A is not a nonpositive integer. The following subprograms are called 
by GAMMF: GSERES, GCHEB, GFRAC, GAMNEG. 


Method : The method was originated by the AEG Computing and Applied Mathematics Center, Courant 
Institute of Mathematical Sciences, New York University. 

(a) If A = 0, 


F(0,X) = Ej(X) = - 


i' + 


log(X) + 

„ n.n.' 
n =0 


(b) If A = -N, for some positive integer N, 


F(-N,X) = ^ 


N 


n; 


Ei(X) - e 


-X Y 


j =0 


X 


(c) If X = 0, 

r(A,0) = r dp = F(A) 

which is the complete gamma function. 
A rational Chebyshev approximation is used: 


(d) For A =it 0, X < ^A + 1, 

• 00 

F(A,X) = F(A) - X^ ) - 

(A + n)n.' 


n =0 


(e) For A:^0, XS\^|a1+'1, 


r(A,X) = e’^X^ 



1 - A 1 2 - A 2 
1 + X + 1 + X+ ■ ■ 
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Accuracy : Complete gamma function: 

Test 1: A = 0.1(0. 1)0.9 by formulas as in reference (a) 

Test 2: A = 1.1(0. 1)1.9 and 10.0(10.0)110.0 

Incomplete gamma function: 

Test 1: A = 1.0(0. 1)2.0, X = 0.1(0. 1)0.9 by formulas in reference (a) 

All test results as compared with table entries of reference (a) were good to about 10 decimal places. 

Reference : (a) Davis, Philip J.: Gamma Function and Related Functions. Handbook of Mathematical 
Functions, Milton Abramowitz and Irene A. Stegun, eds., Nat. Bur. Stand. Appl. Math. Ser. 55, U.S. Dep. 
Com., June 1964, pp. 253-293. 

Storage; GAMMF 610ft locations. 

' > 

Coding information : GAMMF itself is a branching function which according to the values of A and X 
calls the following functions: 

(a) GSERIES(A,X), which computes 

y 

(A + n)nl 
n=0 

(b) GCHEB(A), which computes by a rational Chebyshev approximation F(A) 

(c) GFRAC(A,X), which computes the continued function for F(A,X) 

(d) GAMNEG(IA,X), which computes F(A,X) when A is a negative integer lA 

(Because of the representation Of numbers in the CDC 6600, of A = -N ± e, where e > l.E - io, then A 
is taken to be a negative integer.) 

Subprograms used : System library functions EXP, ALOG. 

Function date: August 1, 1968. 
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Subroutine ITR2 


Language : FORTRAN 

Purpose : Given F(X) = 0, to find a value for X within a given relative error, epsilon, in a given 
interval (a,b). 

Use : CALL ITR2(X,A,B,DELTX,FOFX,El,E2,MAXI,ICODE), where 
X The root 

A The lower bound on X; this value is used by ITR2 as an initial guess 

The upper bound on X; this value is used by ITR2 as a final guess if the entire interval is . 
scanned 

AX, the size of the scanning interval 
The name of a function subprogram to evaluate F(X) 

Relative error criterion 
Absolute error criterion 

A maximum iteration count supplied by the user 

An integer supplied by ITR2 as an error code; this code should be tested by the user on 
return to the calling program: 

0 normal return 

1 maximum iterations are exceeded 

2 DELTX = 0 or negative 

3 a root cannot be found within the given bounds 

4 A > B 

Restrictions : Make A < B, AX positive. A function subprogram with a single argument X must be 
written by the user to evaluate F(X). The name of this subprogram, FOFX, must appear in an 
EXTERNAL statement of the calling program. 

Method : The given function F(X) is evaluated at a given starting point a and at intervals of a specified 
AX thereafter, up to and including a specified end point b. A change of sign of the function across a 
AX interval indicates a possible root in that interval. The interval is then halved successively toward 
F(X) = 0 until the prescribed accuracy is satisfied. The given function F(X) is evaluated once for 
each halving step. 

If the given function is expected to have more than one root between the prescribed starting and end 
points, if is suggested that a sufficiently small value of AX be given so that no more than one root is 
present within a AX interval. A normal return is given upon the location of the first root from the 
starting point a. Additional roots must be located by new entries into the subroutine using a new start- 
ing point a which is just beyond the previous root. 


B 

DELTX 

FOFX 

El 

E2 

MAXI 

ICODE 
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Accuracy : The iteration process is continued until either of two convergence criteria is satisfied, 
criteria are 


It ]X,|><„ 


5 e, 


If XiSei, 




Storage : 260g locations. 
Subroutine date: August 1, 1968. 


These 
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Subroutine OPENMS 


Language : COMPASS 

Purpose : To open a random- access file. 

. Use : CALL OPENMS(U,IX,L,P), where 
U The logical unit number 

IX The first word address of the index 
L The length of the index 

P 0 for numbered indexing; 1 for named indexing 

Restrictions : OPENMS must be the first operation on a random-access file. The file must be a disk file. 
For n index entries, the length of the index must be at least 2n + 1 if using named indexing, whereas 
the index length must be at least n + 1 for numbered indexing. 

Method : OPENMS sets the first word in the index to a positive number for numbered indexing or to a neg- 
ative number for named indexing. The random-access bit, index address, and index length are set by 
OPENMS into the FET of the file for system communication. If the file already exists, the master index 
is read into central memory. 

Storage : lOSg locations. 

Subprograms used : System library subprograms GETBA, SIO$, SYSTEM. 

Error messages : (1) UNASSIGNED MEDIUM FILE XXXXXX 

(2) FILE DOES NOT RESIDE ON A RANDOM ACCESS DEVICE, XXXXXX 

(3) INDEX BUFFER IS OF INSUFFICIENT LENGTH XXXXXX 

XXXXXX is the file name. Termination is abnormal in each case. 

Subroutine date: March 29, 1971. 
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Subroutine READMS 

Language : COMPASS 

Purpose : To read a record on a random-access file! 

Use : CALL READMS(U,FWA,N, I), where 
U The logi car unit number 

FWA The central memory address of the first word of the record 

N The number of words of the record to be transferred 

I . The record number or record name depending upon the indexing mode set by the initial call 
to OPENMS 

Restrictions : The file'must have been opened by a call to OPENMS. 

Method : The disk address of the record is determined using the index. If n words are requested to be 
transferred and there are m words in the record, where m S n, m words are transferred. If 
m > n, n words are transferred. 

Storage : 13 Ig locations. ... 

Subprograms used : System library subprograms GETBA, SYSTEM, SIO$. 

UNASSIGNED MEDIUM FILE XXXXXXX 

FILE WAS NOT OPENED BY A CALL TO SUBROUTINE OPENMS 
RECORD NAME REFERRED TO IN CALL IS NOT IN THE FILE INDEX 
■READ PARITY ERROR* 

SPECIFIED INDEX IN THIS MASS STORAGE CALL .GT. MASTER INDEX OR IS ZERO 
Termination is abnormal. 

Subroutine date : March 29, 1971. 


Error messages : (1) 
(2) 

(3) 

(4) 

(5) 
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Subroutine RECIN 


Language: COMPASS 


Purpose : To read binary records written by the subroutine RECOUT(Jl. 1). 


Use: 1. Type 1 — 

Individual elements (not arrays): 

CALL RECIN(LUN,IT,ICOUNT,Ll,L2,. . .LN), where 

LUN 

Logical unit number 

IT 

Type, equal to 1 

ICOUNT 

Location reserved by the user; RECIN will store the following information in this 
location: 0, end-of-file; nonzero, number of words actually in the logical record; 
if the end-of-file flag was written by a call to RECOUT with lEOF = 1, then end- 
of-file testing must be done by testing ICOUNT for 0; if the end-df-file was written 
by an END FILE statement, then testing for end-of-file must be done by the 
IF(EOF,LUN) statement 

LI,L2,. . .LN 

Individual list elements 

2. Type 2 - 

Arrays: 

CALL RECIN(LUN,IT,ICOUNT,ARRAY,IFIRST,ILAST,INC), where 

LUN 

Logical unit number 

IT 

Type, equal to 2 

ICOUNT 

0, end-of-file; nonzero, number of words actually in the logical record (See ICOUNT 
under type 1) 

ARRAY 

Array name 

I FIRST 

First subscript 

ILAST 

Last subscript 

INC 

Increment 


Examples : 1. CALL RECIN(1,1,K,A,B,ARRAY(1),ARRAY(2)). 

Read a record from logical unit 1 into A, B, ARRAY(l), and ARRAY(2). Note that if the 
record contained only three words, K would equal 3 and ARRAY (2) would be unaltered. 

2. CALL RECIN(1,2,K,ARRAY,1,39,2). 

Read 20 words from logical unit 1 into ARRAY(l), ARRAY(3), . . ARRAY(39). 

Restrictions : If RECIN is used on a file, the only other FORTRAN statements which may be used on that 
file are REWIND and IF(EOF,i). 

The buffer size must be at least 200 13. 

RECIN must be used to read files written by RECOUT and only by RECOUT. 
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Method : RECIN reads into a central memory buffer physical records written by RECOUT, then passes to 
the user the requested logical record via a list giving the elements of the desired logical record. RECIN 
is analogous to a FORTRAN binary READ statement. 

Storage : 30 Ig locations. 

Other coding information : Day file diagnostics and their meaning: 

1. UNASSIGNED FILE MEDIUM FILE TAPEnn - No FET exists for this file. 

Every file has a file environment table that contains information 
describing the file to the system. This error would probably result 
because the file was not defined in the PROGRAM card or the user 
accidentally overwrote portions of his program. 

2. BAD TYPE — The IT parameter was not 1 or 2. 

3. UNCHECKED END FILE — The program attempted to read past EOF with- 

out testing for EOF. 

4. READ/WRITE SEQUENCE ERROR — An attempt was made to read after 

writing. 

Subroutine date: September 22, 1968. 
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Subroutine RECOUT 

Language : COMPASS 

Pbirpose : To write short binary records on a disk or tape in an optimum manner to increase peripheral 
processor and central processor efficiency. These records are to be read by RECIN(I1.1). 

Use : RECOUT may be used for either tape or disk files. 

1. Type 1 - Individual elements (not arrays): 

CALL RECOUT(LUN,IT,IEOF,Ll,L2,. . ,,LN), where 

LUN Logical unit number 

IT Type, equal to I 

lEOF Equal to 1 if an end-of-file flag is desired, otherwise it must be zero. There are two 

methods by which the user may end his file. One method is to call RECOUT with 
lEOF = 1 when the last data record is written. This will cause an end-of-file flag 
(a short length record of less than 512 jq CM words) to be written. RECIN is pro- 
gramed to sense this and will set ICOUNT = 0 when sensed. If this method is 
used, the user must set lEOF = 1 when outputting his last data record since 
RECOUT should not be called with an empty list. For all other calls to RECOUT, 
lEOF must be set to 0. The other method of ending the file is to use the END 
FILE statement. This is the most convenient way of ending the file 
• 

L1,L2,. . .LN Individual list elements 

2. Type 2 - Arrays: 

CALL RECOUT(LUN,IT,IEOF,ARRAY,IFIRST,ILAST,INC), where 


LUN 

Logical unit number 

IT 

Type, equal to 2 

lEOF 

Equal to 1 if an end-i 
under type 1) 

ARRAY 

Array name 

I FIRST 

First subscript 

I LAST 

Last subscript 

INC 

Increment 


Examples : 1. CALL RECOUT(1,1,0,A,B,ARRAY(1),ARRAY(2)). 

Write a record on logical unit 1 containing A, B, ARRAY(l), ARRAY(2). 

2. CALL RECOUT(1,2,0,ARRAY, 1,20,1). 

Write a record containing ARRAY(l) through ARRAY(20). This is equivalent 
to WRlTE(l) (ARRAY(I), I = 1, 20). 
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Restrictions : If RECOUT is used on a file, the only other FORTRAN statements which may be used on that 
file are REWIND and END FILE. 

The buffer size must be at least 20OI3. A normal FORTRAN buffer is this size. 

Files written with RECOUT must be read with RECIN. 

If the list to be written in a logical record is larger than Slljo CM words, then RECOUT offers no 
advantage and should not be used. 

If the programer wishes to write a file containing multifiles using RECOUT, then he must end each file 
by setting lEOF = 1 and not by using the END FILE statement. Consequently, he should then test for 
end-of-file in RECIN by testing ICOUNT for zero. 

Method : Under the CDC SCOPE 3.0 operating system, each binary write commanded by the FORTRAN 
statement WRITE (LUN). . . causes one or more physical records to be output to either a disk or tape 
file. If the logical record size written by the programer is small and the number of records processed 
is large, then excessive usage of I/O routines and equipment results. To decrease this I/O time, 
RECOUT blocks binary data into an optimum record size (512 ]^q CM words) in a central memory buffer 
before transmitting it to the actual disk or tape file. ■ 

Storage : 320 b locations. 

Other coding information : Day file diagnostics and their meaning: 

1. UNASSIGNED FILE MEDIUM FILE TAPEnn - No FET exists for the file. 

Every file has a file environment table that contains information 
describing the file to the system. This error would probably result 
because the file was not defined in the PROGRAM card or the user 
accidentally overwrote portions of his program. 

2. BAD TYPE - The IT parameter was not 1 or 2. 

3. BUFFER TOO SMALL — The buffer size was less than 200l8- 

4. BAD PARAM COUNT — The number of parameters in the call was illegal. 

5. WRITE/READ SEQUENCE ERROR - A write request was made after a 

read request. 

Source : CDC. 

Subroutine date : September 23, 1968. 
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Subroutine SIMEQ 

Language : FORTRAN 

Purpose : SIMEQ solves the matrix equation AX = B where A is a square coefficient matrix and B is a 
matrix of constant vectors. The solution to a set of simultaneous equations and the determinant may be 
obtained. If the user wants the determinant only, use DETEV for savings in time and storage. 

Use: CALL SIMEQ (A, N, B, M, DETERM, IPIVOT, NMAX, ISCALE) 

A A two-dimensional array of the coefficients. 

N The order of A; 1 £ N g NMAX. 

B A two-dimensional array of the constant vectors B. On return to calling program, X is 

stored in B. 

M The number of column vectors in B. 

DETERM Gives the value of the determinant by the following formula: 

DET(A) =(l0l00)^®^^^®(DETERM) 

IPIVOT A one -dimensional array of temporary storage used by the routine. 

NMAX The maximum order of A as stated in dimension statement of calling program. 

ISCALE A scale factor computed by subroutine to keep results of computation within the floating- 
point word size of the computer. 

Restrictions : Arrays A, B, and IPIVOT are dimensioned with variable dimensions in the subroutine. The 
maximum size of these arrays must be specified in a DIMENSION statement of the calling program as: 

A (NMAX, NMAX), B (NMAX, M), IPIVOT (NMAX). The original matrices, A and B, are destroyed. 
They must be saved by the user if there is further need for them. The determinant is set to zero for 
a singular matrix. 

Method : Jordan's method is used through a succession of elementary transformations: Jn-li • • •> ^1- 

If these transformations are applied to a matrix B of constant vectors, the result is X where AX = B. 
Each transformation is selected so that the largest element is used in the pivotal position. 

Accuracy: Total pivotal strategy is used to minimize the rounding errors; however, the accuracy of the 
final results depends upon how well -conditioned the original matrix is. 

Reference : (a) Fox, L.: An Introduction to Numerical Linear Algebra. Oxford Univ. Press, c. 1965. 

Storage : 432g locations. 

Subroutine date: August 1, 1968. 
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Subroutine WRITMS 

Language : COMPASS 

Purpose: To write a record on a random-access file. 

Use : CALL WRITMS(U,FWA,N,I), where 

U Logical unit number 

FWA Central niemory address of the first word of the record 

N Number of central memory words to be transferred 

I Record number or record name depending upon the indexing mode set by the initial call to 

OPENMS 

Restrictions: The filu must have been opened by a call to OPENMS. 

Method: The specified record is written on the file and an address entered in the index to reference the 
record. 

Storage : 102g locations. 

Subprograms used : System library subprograms GETBA, SYSTEM, SIO$. 

Error messages : (1) UNASSIGNED MEDIUM FILE XXXXXXX 

(2) FILE WAS NOT OPENED BY A CALL TO SUBROUTINE OPENMS 

(3) INDEX BUFFER IS OF INSUFFICIENT LENGTH 

Subroutine date: March 29, 1971. 
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SOURCE USTING OF PATS 


OVERLAY( PATS.O ,01 

program MAIN( in put, output ,TAPEl,TAPti>=I N PUT , T A P E 6 = 0UT P UT , T AP E 7 , 
1TAPE8, TAPE9 ) 

COMMON CHAIN (A 102 I 

CHAIN MUST BE OIMENSIONEi) A*NHAX<-6 

KNDEX MUST BE DIMENSIDNEO NBCMAX ♦■2 
DIMENSION INDEX(129) ,<NDEX(802l 

COMMON/BLK 1 / S T A RT T , I T E MT , NB LK , I POw2,NCH , NPRINT, IPLOTA , I PLUTO , OLE SC 
1 AL, DELTA T, SN, NR SKIP, LAP ,NC ROSS, I C K JSS I 2 , 20 » , NC H P , Y L AbE L ( ) , I k,lNUQ,< 

I ,F1 ,F2, I TYPESP , ^^CON,^iE SKIP , IEE,LAbl ,L/.04; 

C0MM0N/61K2/ICH ( 1 A I , CHSUM { 14 ) , N OFF ( 14 ) , CHSUMSDI 14) , SI GMp ( 14) ,RMS( 1 
14) , MEAN( 14) ,SC ALFAlI 14) , C HSUMl ( 1 4 ) , TkAC K ( 1 4 ) , ICHANI 14) ,IFlLT£k(14) 
COMMON/ BLK3/NPT ,TMAX .NPT02, NSPCT , UELF , N64 , NP T 0 1 2 d 
1 , INZERO, MPEAn.NPTliT 
C0MM0N/8LK4/NMAX,MCFiMAX ,NBCMAX 

C0MM0N/BLK5/PCTC .NBINS , UMAX! 1 4 ) , DM I N ( 14 ) , OB i N ( 1 4 ) , B IN S ( 1 DO , 14 ) 
1,CH1S0C 

C0MM0N/BLK6/ISAVEt)4, IR I , I XPLOT , 1 0 ATA , I 2 , I SPECT 
C0MM0N/BLK7/I w'ORC( 11 ) 

COMMON/ BLK8/ I AUTOSPI 14) , I AUTOCOI 14),ICRSPI20),I CRC0R(2 0) , 1TKA(20» , 
1ICOH(20) 

CGMMOF7BLK9/NFILTP,FREOF(SO),WGHTF(SC ) 

C OM MON/ B LK 10/NR CP07 
DATA NRCRD7/0/ 


NMAX=1C24 

NCHMAX=14 

NBCMAX=80'" 

IZ=IRI=1 

IXPLOT=IRI 4NMAX44 
I SPECT = I SA VE64= IXPL0T4NMAX42 
IDATA=IS AVE64464*NCHMAX 
PATS=4HP ATS 

CALL OPEMMS (9, INDEX, 129,0 ) 

CALL OPENMSId, KNDEX, NflCMAX42,0) 

PRINT 1000 

,1000 FORMAT! IHl //////////23X 1H*43( 2H ♦ ) /2 3x1 H*85X 1H«/20X IH-l-ax^P R 0 G R 
1AM FOR ANALYSIS OF TIME SERIE S^dKiH* 
2/23X IH’FSSXIH*/ 23XlHt'SX*C .G. BROWN , T.J. BROWN, AND J.C. HARDIN FOP A 
3C0USTICS DIVISION, NASA-LRC, 197 3-3X 1 H* / 23X lH‘i'd 5X 1 H* /23XlH*43i2H 

4A) ) 

READ NAMELIST AND CARD INPUT 
100 CALL OVE RL AYIPATS, 1,0) 

READ ONE BLOCK CE INPUT TIME SERItS DATA FOR NCH CHANNELS FkOM 
BINARY TAPE AND STORE ON RANDOM ACCESS FILE 
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COMPUTE FOURIEK TRiViMSPORM FOR NCri CHANNELS ANU STORE ON RANDOM 
ACCESS F !L£ 

REPEAT UNTIL ALL DATA IS READ AND PROCESSED. , 

CALL OVERLAYIPATS ,2,0) 

AVERAGE NBLK POWER SPECTRA FOR EACH CHANNEL FOR AUTO SPECTRA AND 
APPLY spectral FILTERo COMPUTE AUTUCOR3 ELAT I ON. PRINT FANFOLD PLOTS 

CALL OVEHL AY(PATS,3,0) 

TO COMPUTE CROSS SPECTRA 

GET NBLK TRANSFORMS FOR BOTH OF THE CHANNELS IN EACH SELECTED PAIR 
AVERAGE THE PRODUCTS FUP NARROW BAND SPECTRA. COMPUTE C KOS SC OkRc LA T i QNS , 
COHERENCE, AND TRANSFER FUNCT IONS . PR INT FANFOLD PLOTS. 

IF(NCRGSS.GT.O) call OVERLAYI PATS,'*,0) 

GO TO 100 
END 

SUBROUTINE' PL OT NB ( YL A B E L , FR AM EL , NF , X PLU T , k ! , .i SP C T , I LOb E i , F 2 , PL A BE 
IL.NP, IFF, T SEARCH) 

DIMENSION YLAriEL(2),FKAMEL(l),XFLOTU),RIll),PLAbuL(SI,l-ElD4S), 

1IDEN(6) 

ILOG - CODE FOR TYPE OF SCALE FOR AXES 
= 0 BOTH scales LINEAR 

= 1 HORIZ. SCALE LOG, VFRT, SCALE LINEAR 
■= 2 HURIZ. SCALE LINEAR, VERT. SCALE LOG 
= .3 BOTH scales LOG 

OPTIilMS I AND 3 ARE NOT USED BY PATS 

ILQGP1=T LOG+l 

IF( I SEARCH, EO, 0 ) GO TO 6 

Jl=l i J2=NSPCT 

DO 1 I = l ,NSPCT 

J=I 

IF( XPLOT ( I ) .GE.Fl) GO TO 2 

1 CONTINUE 

2 Ji=J i IF( ILOG.ED.l.OR. ILOG.b J.3) J 1 =MAXC ( J I , 2 ) 

00 3 I=J1,NSPCT 

J = NSPCT-H-J1 

IFI XPLOT (J ) .LF.F2) GO TO 4 

3 CONTINUE 

4 J2=J 

NFi\ = NF/10 $ IF(NFw*IO,LT.NF) NFw = NFw*-i 
NFP=NF/10 t TF (NFPk10.lt. NP) NFp=NFP+i 
IF ( J1.LT.J2) GO TO 3 

PRINT 900, (PL ABEL( I ) , I =1 ,NFP ) , ( FRAMLLI I ) , 1 = 1 ,NFW) 

900 FORMAT ( K-OBANDW I CTH FOR PLOTS TOO NARRuw*/k NO PLUi' GE,NERAT£D FOk * 
llOAIO) 
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RFTURN 

5 'NPL'JT=J2-J H-1 

lF(NPLOToC e.NiSPCT I GO TQ ICC 
C MOVE plotting REGION! TO BEGINNI'IG OF ARRAY 

00 4002 ! = J 1, J2 
J=I-JU-1 

4002 XPLOT( J ) =XPLOT ( n 

6 CONTINUE 

DO 4003 I=J1,J2 

J=!-JH-1 

RI ( J )=RI ( I) 

4003 CONTINUE 
100 CONTINUE 

103 GO TO ( ICE ,105 , 106,106 I , ILOGPl 
C LINEAR VERTICAL SCALE 

1C5 CALL ASCALEIRI , 10,,NPL0T,1,10.) 

YMIN=RI (NFLOT 4-1 ( % Y M AX = Y M I N 10. «R H NPL OT f 2 ) 

GO TO 107 

C LOG VERTICAL SCALE 

106 SMAX=-10C. 

on 103 I=1,NPL0T t n-IRKII) 112,112,113 

112 RT( I )=-iOO. i GC TO 103 

113 RI( I ) =ALPG10( K I ( n ) $ SMAX = AMAX1 ( PI ( I ) , SNAX ) 

108 CONTINUE 

IMAX=SMAX 

IE ( IMAX, LIo SMA X ) IMAX=INAX«-1 
IMIN=IMAX-5 

YMIN=RI (NPL0T4-1 » = IMiri $ YMAX=IMAX $ R I ( NP LOT 4- 2 ) = . 5 

00 114 I=l,NPLOT 

IF(RI( I I.LT.YMIM RI(n=YMIN 

114 CONTINUE 

107 CONTINUE 

IF( IFE.EO.O I GU TO 24 
IF(NPL0T.GT.256) NPLQT=256 

FFUm) = YLA6tL(l) $ FF I 0 ( 2 ) = YL A3 EL ( 2 I $ F F 1 0 ( 3 ) =FR A ME L ( 1 ) 

FFID(4) = FRAMEL ( 2) t FF I 0 ( 5 ) =FR A ME L( 3 ) 

IF( ILOG.GT. 1 I ENCOOE ( 54, vOl , IDEN I PLABEL 
IF( ILOG.LF. U FNCOOE (60,9C2, lOENJ PLABEL 
902 FOKM ATI 5X , 5A1C , EX ) 

CALL HANFOLUIR I ,NPLP'^ , i , 1 , NPLOT, lUEN , Ih * , 1, , YM A X , YiM I N , F F 10 , NF W4- 2 , 
U32,0,XPLOT ,9hFREGUENCY 1 
. 901 FORMAT! >f^LOG *5A10) 

24 RETURN 
ENO 

sue ROUT INF FANEOLO(PLOT,NPT,K,,NF,NMAX,iuEN,CHAR,PNOKM,PMAX, PM IN, 
lYLABFL ,NYL , LINE , I WRI T E X , X ARR A Y , X L AtJE L i 
DIMENSION PLOT (NMAX, 1) , C H AR I 1 I , I OEN ( 5 , 1 ) ,Pi40RK( 1) ,YLABELl2) 
1,PMAX(1 ),PMIM1 J ,PSCALE I 10) , PL INEI lit I ,XAkRAY I 1 ) 

IF(LINE-12C) 7,7,0 

7 NCHAR=104 f GO TO 9 
a NCHAR=124 
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9 CONTINUE 

IF(NF.GT.IO) NF=1C 

PRINT 900, (YL/idEL( I) , I = 1,NYLI 

900 FORMM ( * l*43X8a 10) 

IF( IWRITFX.EQoO ) GO TO 10 

PRINT 907, (XL ABEL , J = i, 0 ) , ( J, X4RR aY( J ) , J = 1 ,NPT ,K ) 

9C7 FORMAT ( //^IXoThE *A10,* CLiUE SCALE IS AS FUL L OW S*/ 5X «C OUh KAlC.'tl 
16X*CGUE *A10)/(IB,E1R. b,IB,E19o5,Io,cl4.b,Id,E19. S,lB,bl9.5)) 

1C CONTINUE 
NSKIP=0 

IF( NPT/K-15 ) 16,16,17 

16 NSKIP=2 $ GU TO 13 

17 IF(NPT/K-25I 18,16,13 

18 NSKIP=1 
13 CONTINUE 

DO 1 1=1, NF 

IFIPMINt I loNEoPMAXd ) ) GO TO 1 
PMIN( I )=PMAX( I ) =PLOT (1,1) 

Dll 2 J = 1 ,NPT,K 

PMINd ) = AMIN1(FMIN(I ),PLOT(J,D) 

PMAX( I ) = AMAXl ( PMAX( I ) , PLOT (J, 1 ) ) 

2 CONTINUE 

1 PSC ALE ( I ) = ( FMAX ( I )-PM IN ( I ) )«PNOTM( I ) / (NCHAR-',. ) 

PRINT 90 1, ( I, CHAP ( I ) , ( IDEN( J, I ) , J = l, b ) , PNORf-K I ) , PM AX ( I ) , pM IN ( I ) , 
IPSCALEl I ) , I =1 ,NF ) 

901 FORMAT ( ///Sax^PLOT OESCR I PT I0 n«/1' FUNCT 1 0W*6 7X*SC A L E* / ■!“ ND.<‘4X*C 

IHAR ACTE R *6X*I UEMT IFI CAT ! ON =i‘3 7 X* F ACT OR ■i'9 X ^NAX I MUN’«‘8 X^M I N 1 nU M«6X* RE S 
2CLUT ll.!M«/( I 5,9X ,Al,bX,5A lC,9£ib.8)) 

PLINE(2)=lFi( t, PL INE (NCHAR ) = 1H) $PLlNt(l) = lH 


IE( NCHAR.E0.129 ) PRINT 9C2 
IF(NCHAR,FQ,109 1 PRINT 905 

902 FORMAT! //20X*j. 12 2 3 

13995566 7 788991 

20 10 11 11 12*/ 10X='0. . ..5.. ..C. . . .5. . . .0. . . . 5. . . .0. i . .5. . . 

30Oa««o5oO',oOao««5«oooV^«oo«5aoa«Oo«oa5o«ooO*»«»5».«*O««»*3«*»*O**»* 

•t5....:....5....0*) 

905 F0rMAT( //2Gx*i 1223 

139 9 55 66 7 7(1 8 991 

20*/ 10X*0. . . .5. . ..C.. . .‘3. . . .0. .. .5.. . .0. . . .5. . . .C . . ..5. . ..C. . ..5. . . 
3.0....5....0....5....0....5...<,0..o<,5<>...0*) 


IE(NSK!P.ECJ.O) GP TO 20 
DO 21 1 = 1, NSK ! P 
21 PRINT 903 
906 FORMAT (IH ) 

20 CONTINUE 

OP' 3 J= 1 ,NPT , K 
on 5 I =9, NCHAR 
5 PLINE( 1-1) = 1H. 

OU 9 1 = 1 ,NF 

P=( PLOT( J, I )-PMIi\( I) )/PSCALE( I ) 
IP=P+3. 5 



APPENDIX H 


IF< I P.LE ,2 . .OP . IP.GE.NCH4R > GO TO 4 
IF( PL INt ( » P ).NE .IH » GO TO 6 
PL INE ( I P ) =CHAK (II 
GO TP 4 

6 PLINE ( I P I=1HX 
4 CONTINUE 

PRINT 903, J, ( PLINE( I I , I = 1,NCHARI 

903 FORMAT { I 8 , 124i n 

IF( NSKI P .EO.O I GO TO 3 
un 19 I=1,NSK1P 
19 PRINT 908 
3 CONTINUE 

IF(NCHAR.EQ.124 ) PRINT 904 
IF I NCHAR.ro. 104 I PRI NT 9C6 

904 FORMAT (10X*0.o..5...<.0.,.,8...,0. 


15....0....3....0....b....0....5....0....8....0....3....0....3....0 

3*/ 20a* 1 i 2 2 33 4 

14 5 5 6 6 7 7 8 d -» 910 101 

21 11 12=:‘) 

906 FORMAT! 10X*0....6....0....5....0....6....0. <,..3....0....b>. ...0. ... 

13. . . .0. . . ,5 .. . . 0. . .. 5. . . .0. . . .5. . ..0. .. . b.. . . /20X*1 i 2 

22 3. 3 4 4 b b 6 o 7 7 8 d 

39 9 10*1 

RETURN 
END 

SUBROUTINE FOUR T ( OAT A , NN , NO I M , I S I GN, I FORM , WORK ) OuO^ 

C 00w3 

C the COOLEv-TUKEY fast FOURIER TRANSFORM IN USASl BASIC FURTR ' 0uC4 

C 0 C C b 

C TRANSF0RM(K1,K2 = SUM ( DA T A ( J 1 , J 2 , . . . I* E XP ( I S I GN*2 *P I *SQ 00C6 

C *( ( J 1-1 ) *(K 1-1 ) /NN(1 1 «■( J2- 1 l*(K2-l)/N.N(2 I »•. .. I) ) , SUMMED FOR 0LJ7 

C Jl, K1 BETWEEN 1 AND NNdl, J2, K2 BETWEEN 1 AND NN(2), ETC. 0008 

C THERE IS NO' LIMIT TO THE NUMBER Of SUBSCRIPTS. UATA IS A 0CC9 

C MULTIDIMENSIONAL COMPLEX ARRAY WHOSE REAL AND IMAGINARY OOiO 

C PARTS ARE ADJACENT IN STORAGE, SUCH AS FORTRAN IV PLACES THt „011 

C IF ALL IMAGINARY PARTS ARE ZERO (UATA ARE DISGUISED REAL), S OCl^ 

C IFORM TO ZERO TO CUT THE RUNNING TIME BY UP TO FORTY PLkCENT i,C13 

C OTHERWISE, IFCJRM = 1 . THE LENGTHS OF ALL DIMENSIONS ARE 

C STORED IN ARRAY NN, OF LENGTH NDIM. THEY KAY BE ANY POSIT IV OCib 

C INTEGERS, THO TFE PROGRAM RUNS FASTER ON COMPOSITE INTEGERS, OGlo 

C' ESPECIALLY FAST ON NUMBERS RICH IN FACTORS OF TWO. ISlGi'i IS 0017 

C OR -1. IF A -I transform IS FOLLOWED BY A «■ 1 ONE (OR A H CC18 

C BY A -1) the original DATA REAPPEAR, MULTIPLIED BY NTOI (=NN 0019 

C NN(2)*...I. TRANSFORM VALUES ARE ALWAYS COMPLEX, AND AkE RE 0020 

C IN ARRAY DATA, REPLACING THE INPUT, IN ADDITION, IF ALL CG21 

C DIMENSIONS ARE NOT POWERS OF TwO, ARRAY WORK MUST BE SUPPLIE 00^2 

C COMPLEX OF LENGTH EQUAL TO THE LARGEST NON 2**K. DIMENSION. C023 

C OTHERWISE, REPLACE WORK BY ZERO IN THE CALLING SEQUENCE. 0C24 

C NORMAL FORTRAN DATA ORDERING IS EXPECTED, FIRST SUBSCRIPT VA cO^-b 

C FASTEST, ALL SUBSCRIPTS BEGIN AT ONE. d026 

C uC2Z 
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C RUNNING TINE IS MUCH SHORTER THAN THE NAIVE NTJT**2, bElNG C02rf 

C . GIVEN BY THE EOLLOWING FORMULA. OECUMPOSE NTDT INTO 002V 

C 2«*K2 * 3**K3 * 5«*K5 * .... LET SUM2 = 2*K2, SUMF = + 0C30 

C ... AND NE = K3 KS .... THE TIME TAKEN BY A MULTI- OEol 

C DIMENSIONAL TbANSEOkM ON THESE NTOT DATA I S T = TO h NT0T*IT 0032 

C T2*SUM2hT 3*SUMF i-TA^NF ) . CN THE CUC 33C u I FL OAT I NG POINT ADD 0G33 

C OF SIX M ICRQSECCNDS) , T = 3000 NTGT ♦ ( 500*- ^ 3<= S UM 2«- 66* SUMF t 003't 

C 320*NFI MICROSECUNOS ON COMPLEX DATA. IN ADDITION, Tht 0035 

C ACCURACY IS GREATLY IMPROVED, AS THE FMS RELATIVE ERRUk IS DC3o 

C BOUNDED BY 3*2 * * ( -B ) * SUM ( F ACT OR ( J I ** 1 . 5 ) , WHERE B IS THE HUM 0C37 

C OF BITS TN the FLCATING POINT FRACTION AND FrtCTORIJI ARE THc 0036 

C PRIME FACTORS OF NTOT. OCjv 

C U040 

C THE DISCRETE FOURIER TRANSFORM PLACES THREE RESTRICTIONS uPO 0046 

C DATA. 0049 

C I, THE NUMBER OF INPUT DATA AND THE iviJMbER OF TRANSFURM VAL 005u 

C MUST BE THE SAME. 0051 

C 2. BOTH the input DATA ANu THE TRANSFORM VALUES MUST KEPRES 005,: 

C EQUISPACED POINTS IN THEIR RESPECTIVE DOMAINS OF TIME AND uC53 

C FREQUENCY. CALLING THESE SPACINGS DELTAT AND DELTAF, IT MUS 0054 

C ■ TRUE THAT DEL T A E = 2*P I / ( NN I I) *DE L T AT ) . uF COURSE, DELTAT NEE 0055 

C BE the same FOR EVERY DIMENSION. 005o 

C 3, CONCEPTUALLY AT LEAST, THE INPUT DATA AND THE TRANSFORM 0C57 

C- REPRESENT SINGLE CYCLES OF PERIODIC FUNCTIONS. 0056 

C 0059 

C example 1. THRFE-DTMEMSIONAL EOkWARD FOURIER TRANSFORM OF A uOoO 

C COMPLEX ABRAY DIMENSIONED 32 BY 25 bY 13 IN FORTRAN IV. OCoi 

C DIMENSION DATA( 32,25, 13) ,WUKK( 50> ,NN(3) 0C62 

C COMPLEX DATA 0063 

C DATA NN/32,25,13/ ' 0064 

C DO 1 1=1,32 0C63 

C on 1 J = l ,25 . 0066 

C DP 1 K=l,13 0067 

C 1 DATA! I , J,K) =COMFLLX VALUE uCo6 

C call F0URT(DATA,NN,3,-1, 1,W0Rk) G0C9 

C 0070 

r EXAMPLE 2. QME-U!MLNSI0NAL FURWARD TRANSFORM OF A REAL ARRA 0C71 

C length 64 IN lORTRAM II. 0072 

C DIMENSION DATA(2,64) u073 

C IDO 2 1 = 1,64 0074 

C DATA( 1 , I l=REAL PA^T 0075 

C 2 0ATAI2, I ) = C. u0 7o 

C^ CALL F0U0T(nATA,64,l ,-l ,0,0) 0077 

C- 0076 

DIMENSION DAT A( 2) ,NN I 1 ) , II ACT( 32) ,wDRK( 1 ) 0079 

Wl=1.00 0C6u 

wR=1.00 0081 

9 ST PR = 1. CO oCb2 

RSTPI=1.C0 OOtt-. 

TW0PT=5.2B3185307 0C64 

IF(NDIM-1)920, 1,1 o065 
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1 NT0T=2 QCoo 

00 2 QCtil 

lF(NNt IDIM) )920,92 j,2 OOBo 

2 NTnT = NTOT>^NN( TDIM) 0Co9 

C U09u 

r, MAIN LOUP PriH tACh OIMtivjSlUN OOvl 

C 0C92 

iNPl = 2 0C9P 

on 910 ro!M=l ,NOIK 0099 

N = N,N(IOIH) UC9i3 

NP2 = NPl=i'N CG96 

IF(N-1)920,90C',5 009/ 

C 009o 

C F AC TOP. N 0C99 

f. 0100 

5 M=N Old 

NirwO = NPl 0102 

1F=1 0103 

IDIV=2 01o4 

10 I iJUQT = M/ lUI V OiCb 

IR£M = H-I OT V^UOUCT 0109 

IF( IQUOT-1 OIV ) 5Ctil, 11 0107 

11 IF( IR£M ) 20 , 12 , 2C OiOd 

12 NTW0 = 'NTw3*-NTwn ol09 

M=IQUOT Clio 

GO TO 1C 0111 

20 IUIV=3 uli2 

30 IQUOT = M/Ii)IV 0113 

[REM = M-i,ji v«U:UCT 0114 

IF( IOUnT-rO!VI6C,31, 31 • 0113 

31 IF( IPBM)9(',32t4C 0116 

32 IFACT(IF)=I0IV 0117 

IF=IF + 1 Cl la 

M=Tgu<^T 0119 

GO TO 3C 012C 

4C ID!V=!I)IW2 0121 

GLf TO 30 . ■ 0122 

50 IFl IRFM ) 60, 51 , 6C 0i23 

51 NTWO=MTwO+NTwn 0124 

GO TO 70 0U5 

60 IFACT(1F)=M 0126 

C ■ 0127 

C SEPARATE FOUR CASES-- 0123 

C 1. COMPLEX TRANSFORM (JR REAL TRANSFORM FOR ThL 4Th, :>TH,E 0129 

C DIMENSIONS. 0130 

C 2. real ■transform for the 2ND CR 3RD OIMEnSICN, METHOD — 0131 

C TRANSFORM FALF THE DATA, SUPPLYING THE OTHER HALF BY C 01j2 

C JUGATE SYMMETRY. 0133 

C 3. REAL transform FOR THE 1ST DIMENSION, N ODD. METHOD — 0134 

C TRANSFORM HALF THE DATA AT EACH STAGE, SUPPLYING THE D 0135 

C HAIF BY CCNJUGATE SYMMETRY. 013t> 
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C 4. REAL TRANSFORM FOR THE 1ST DIMENSION, N EVEN. METHOD- CiJ/ 

C TRANSFORM A COMPLEX ARRAY OF LENGTH Iv/^ WHOSE REAL PAR 0136 

C ARE the even numbered REAL VALUES AND WHOSE IMAGINARY 0139 

C ARE THE 000 NUMBERED REAL VALUES. SEPARATE AMU Si^PPLY 01‘+0 

C the SECCND half by conjugate symmetry. 0141 

C 0142 

70 NUN2=NP1* ( NP2/NTWO I 01^3 

ICASE=1 0144 

IF( IOIM-41 71,9G,9C ol4S 

71 IF( IFORM 172,72, 9C ul4o 

72 ICASE=2 ’ ul47 

I F{ I niM-1 173,73,90 0i96 

73 ICASE=3 0149 

IF (NTWiJ-NPl )9C , 90,74 0130 

74 ICAS£ = 4 • ul-3i 

NTWO = MTwO/2 OlD^; 

N= N / 2 ' o i 3 3 

NP2 = NP2/2 OlS't 

NT0T=NT0T/2 0133 

T=3 0i3o 

DO 8C J=2,NTOT 0157 

DATA(J) = DATA( 1) C136 

80 I = I«-2 0139 

90 I1PNG=NP1 016C 

IF( ICASE-2 ) IOC, 95, ICC OloT 

95 I 1RNG = NP0* ( IRNPPEV/2 ) 01&2 

C U 1 6 3 

C SHUFFLE UM THE FACTORS UE TWO IN N. AS THE SHUFFLING Clo4 

C CAN RE DONE OV SIMPLE INTERCHANGE, MC WORKING ARRAY IS NcEDE Oloa 

C O' 1 6 6 

lOD IF(NTW0-NP1 )6CC,6LD, 110 0167 

110 , NP2HE=NP2/2 0166 

J= 1 01 69 

DO 150 I 2= 1 ,NP2 ,NUiN2 u17l ' 

IFU-12) 120, 130, 130 ' 0171 

120 ! 1MAX=I2+NTN2-2 w172 

DO 125 I li 12, I 1MAX,2 . 017a 

DO 125 I 3=1 1,MTCT,NP2 0x74 

J3-=J4-I3-I2 0x75 

TEMPR=DATA( 13) 0176 

TEMPI=DATA( 13+1 ) 0177 

OATA( I3)=nATA( J3) 0176 

DATA! 13 + 1 ) = DATA (J3H ) Ci79 

•DATA! J3) =TEMPR 0160 

125 DATA(J3+1)=TEMP1 0l6i 

130 M=NP2HF 0162 

140 IF(J-M)150,15C,145 0163 

145 J=J-M 0164 

M=M/2 j165 

IF{ M-M0N2) 150, 140, 140 uloo 

150 J=J+M 0167 
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C Olbo 

C MMN LOOP rUR HACTOF.S OF TWO, pfcRFOR.M FUOkIcR TRANSFORMS OF Olb-i 

C length four, with one of length TwO if NEEUEO. the TwIOL/LE OiSiO 

C W = t XP( I SIGN*2<=PI'^SQRT ( -1 i *M/( 4*MMAX ) ) . CHECK FOR w=IblGHi‘SU 0191 

C AND REPEAT FOR w= 1 S I GN « SQP T ( - 1 ) *CGN J UGA T E ( w J . J19,i 

C Jl9S 

NaN2T=NaN2fN0N2 0194 

I PAR=NTwO/NPl 0199 

31C IF( IPAR-2) 350, 330,320 0196 

32C IPAR=IPAR/4 0197 

GO TO 310 ol9b 

330 DO 340 I1=1,I1RNG,2 0199 

DO 340 J3=I 1,NCN2,NP1 u2u0 

DO 34C K1=J3,NTCT,N0N2T 0^01 

K2=K1+N0N2 0202 

TEMPR=DATA < K2 ) 

TEMPI=DATA (K2H » 020-+ 

DATA (K2 ) =DATA ( K 1 )-TEMPR o2v5 

0ATA(K2+1)=UATA(KU-U-TEMPI 02Ct> 

DATA ( K1 ) =0ATA ( K1 ) HEMPR 029? 

340 0ATA(K1<-11=DATA(K1*-1 1 t-TEMPI 0208 

350 MMAX=N0N2 0209 

360 IF( MMAX-'K|P2HF ) 3 /C ,600, 600 0210 

370 LMAX=MAXO( NHN2T ,MMAX/2 ) 0211 

IF(MMAX-NQN2I 405,405, 380 0212 

380 THETA=-TwnPI*FLrATlNGN2l/FL0AT(4«MMAX) 0^13 

IF( ISIGN>4C0, 390,390 0214 

390 THETA=-THETA 0215 

400 WR=C0S(THETA) 0216 

WI=SIN( THETAI 0217 

WSTPR = -2.’^Wl*Wl 0218 

WSTPI = 2. *WR>l'WT 0219 

405 DO 570 L=N0N2,LMAX,N0N2T 0220 

M=L- 0221 

IF( MMAX-NnN2)42C,42Q,410 , 0222 

410 W2R = wR*^WR-Wl*WI 0223 

W2I=2.*WP’S=W! u224 

W3R=W2R*WR-W2I*WI 0225 

W3I=W2R*WI+W2I *WR 0226 

420 00 530 I1=1,11RNG,2 0227 

DO 530 J3=T 1,M0N2,NP1 02,i8 

KM1N=J3MPAR*N 0229 

IF{MMAX-N0N2)43C,430,443 , 0230 

430 KMIN=J3 C2j1 

440 KDIF= IPAR^MMAX 0232 

450 KSTEP=4«K0IF 0233 

DO 52C K1 = KHTN,M0T,KSTEP u^;34 

K2=KH-K0IF 02j5 

K3=K2<-KD!F 0236 

K4=K3+KDIF 0237 

IFI MMAX-N0N2)46C, 460,480 02j8 
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<^6 0 

U1P. = 0ATA(K1 M-DATA(K2 ) 

0239 


U1 r- DATA («■ l«-l ) ♦■GAT 4 ( K 2+ 1 1 

0240 


U2R = DATA (K3 UGATA(K4) 

0241 


U2I = OATA<K3*l ) ♦CAT-A(K4^1I 

0242 


tJ3R = DATA (K1 I-1 ;atA(K2 I 

0243 


U3I = nATA(KlHt-CATA(:<.2vi> 

0244 


IF( ^SIG^J)47C',47 5,475 

U245 

470 

04P. =OATA (K3Fl ) - CAT a( 1 ) 

0246 


IJ4I =0ATA (K4I-GATA(K3 1 

02 47 


00 TO 510 

0248 

475 

■J4R = GATA (K4M 1 - OAT A( i<3*-l ) 

0249 


U4I=0ATA (K3)-0ATA(K4 ) 

0230 


00 TO 51 C 

0 ii5 1 

480 

T2R = W2R*nAT A( K 2 )-H?I -n) AT A ( K2 1 1 

0252 


T2I=W2R*OATA(K2 + 1)»Fi 2I ’i^UAT A( K2) 

0^53 


T3R-WP *UAT A(K2 ) -W I AT A ( K.3 ♦ 1 1 

0234 


T31=WK*r)AT A (K3 + 1 ) ♦WI 4>DAT A ( K3) 

0253 


T4R = W3R*0ATA( K4 ) ->a3 1 *U A t A ( K4f 1 1 

0256 


T41 = WJR*[)ATA( K4 + 1 )f'w 31 ’*‘UATA( K4) 

0257 


IJIP-OAT A (K 1 1 ♦-T2R 

0^30 


U1I=DATA (K14-1) ♦T2! 

0259 


U2R=T3R<-T4R 

02o0 


02I = T3UT4! 

GzOl 


U3P = r)ATA(Kl)-T2R 

026 2 


L)3I=DATA(^i♦•l»-T2I 

0263 


IF( I SIG.M )490,f'CC,500 

0264 

49 0 

U4R=T31-T4! 

0^65 


U4I=T4R-T3R 

0266 


GO TO 510 

0267 

500 

U4R = T‘tI-T3^ 

0268 


U4T =T3R-T4‘-' 

0269 

510 

DATA(Kl) =iilRfU2P 

02 70 


0ATA(K1^-1)=U11 ♦U2T 

0^7 1 


0ATA(K2 » =U3R*-II4P 

0272 


OATA(K2^-ll=U3I ♦L4’ 

0273 


0ATA(K3) ='ilP-i)2R 

0274 


^DATA(K3fl)^UlI-U2! 

0275 


^)ATA(K4»=U3R-iJ4R 

0276 

5^0 

data (K4<- 1) =U3I-U4I 

6277 


KMiM = Zt«(KMIN-J3)vJ3 

0278 


KOI F = KST EP 

0279 


IF(KniF-RP2)450,b3G,530 

0280 

530 

CONTINUE 

0281 


M=MMAX-M 

0282 


IF( IS!GN)540tR5C,550 

0263 

540 

TEPPP=WO 

0284 


WR=-WI 

02o5 


■WI = -TEMPR 

0<;86 


GO TO 5oC 

0287 

550 

TEMPP=WP, 

0268 


WR = WI 

0289 
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WI=Tt::’^PP 

56C IF( M-L'^AX) 56!)f 565,<^10 

565 TEMPR = '^« 

WP=HR'"AlSTPP,-hH‘kSTP! i-HP. 

570 WI=vVT*WSTPR«-TEMFR*WSTPI ♦■WI 
IPAR=3-IPAR 
MMAX = MMAX«-MMAX 
go' to 360 
C 

r, MAIN LOOP FOR FACTORS NOT fcO'JAL TO TWO. APPLY THE TWIDULE F 

C W=EXP{ 1 SIGN*2*P I*SyRT (-U “M JZ- U JI-J2 ) / (NP2*l EP U» , THEN 

C PERFORM A FOURIER TRANSFORM OF LENGTH IFACTIIF), MAKING USE 

C CONJUGATE SYMMETRIES. 

C 

600 IF(NTWQ-NP2)605,7C0, 700 

605 IFPi=N0N2 
1F=1 

NPlHF=NPl/2 

610 IFP2=IFPl/IFACT(in 
J1RNG=NP2 

IF( ICASE-3 ) 612 ,611,6 12 

611 J1RNG=(NP2MFP1 1/2 
J2STP=NP2/ IFACTI IF) 

J1RG2= ( J2STPtIFP2)/2 

612 J2MIN=H-IFP2 

IF( IFPl-No2)615,6A0, o40 
615 00 635 J2=J2M1N,IFP1 , IFP2 

THETA = -TW0PI*FLCATU2-11 /FLOAT! NP 2 1 
IF( ISIGN1625, 620,620 
620 THETA=-THETA 
625 SINTH=ST NI THETA/2. ) 

WSTPR=-2o*S INTH*S1NTH 
WSTPI=SIN(THFTA| 

WR=WSTPR f 1. 

WI=WSTPI 

J1MIN=J2HFPI 

. 00 o35 J 1 = J IMIN , JlRNCi, IFP 1 
I1MAX=J1 H lRNG-2 
on 6 30 I 1= J 1, I IMAX,2 
DO 63C I 3= I 1 ,NTCT ,NP2 
J3MAX=I 3«-IFP2-NPl 

00 630 J 3= I 3, J 3max,NP1 
TEMPR=DATA (J3) 

DATA! J3 1 =HATA ( JSX-wR -DATA ( J3«-l 1 I 
630 DATA! J3+l) = TcMPR«WH-0ATM Jj<-1)*WR 
TEMPR=WR 

wR=WR4=wSTPR-Wl <‘FSTPI «-WR 
635 WI=TEMPR*WSTPI I-UK-WSTPRVWI 
64 0 THE T A = -T WOP I /FLOAT! [FACT! IF )» 

1 F( I SIGN' I 650,64 5,645 
645 TH£TA=-THETA 


02VO 

0291 

0292 

029 3 

0294 

0295 
0^96 

0297 

0298 

0299 
C300 

0301 

0302 
C3u3 
0304 
o306 
0306 

030 1 

0308 

0309 

0310 

0311 

0312 

0313 

0314 
3315 

0316 

0317 

0318 

0319 

0320 

0321 
0^22 
0323 
03il4 

0325 

0326 
0 32 7 

0328 

0329 

0330 

0331 
03 32 

0333 

0334 

0335 

0336 

0337 

0338 

0339 

0340 
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650 SINTH=SIN( THETA/2.) 0:>^i 

WSTPR = -2.^SIKTH«SIfvlTri 

WSTPI=SIN(THETA) 0A43 

KSTEP=2->^N/ I FACT ( IF ) 0344 

KRANG=KSTEP*( IF ACT ( I F I /2 ) +1 cJ4a 

00 698 Il=i,IlRNG,2 034o 

on 698 I 3=1 It NTCT.NP2 034/ 

DO 690 KMI N=l ,KRA'\iG,i<.STEP 0348 

J1MAX=I34-J iPA'G-lFFl 0349 

00 83C J1=I3, JlMAXtlFPl u35C 

J3MAX=JH-IFP2-NF 1 0-851 

DO 680 J3=Jl, J3PAX,\IP1 0353 

J2MAX = J 3 «■! FPl- I FP2 C35J 

K = KMIf'j4-( J3-JH- ( Jl- 13 I / IFAC T ( IF ) )/v-4PlhF - 0354 

. IF( KMIN-1) 655t655,665 0355 

655 . SOMR=G. Coao 

SUMI=0. ' 0357 

00 660 J2=J3, J2MAX,IFP2 .8358 

SUMR=SUMR«-DATA { J2) 0359 

660 SUM1 = SUM I f OAT A ( J2 + 1) 03o0 

WORKU) = S;j'1R 8 361 

WORK (K.vl ) =SUM1 0362 

GO TO 680 0363 

665 KCQNJ=Kf2*( N-KMIN4-1) 0364 

J2= J2MAX 03o5 

SUMR=DAT A{ J2) 0j6o 

SUMI = 0ATA( J2<-1 ) 036/ 

OL0SR=0. 0368 

0L0SI=0. 0369 

J2=J2-IFP2 037o 

670 T£MPR=SUM9 03/1 

TEMPI=SUMI 0372 

SUMR = TW0WR*SUKR-0LL)SR4DAT A( J2 ) 03/3 

SUMI = TWOWR*SUf«tl-nLOS I tOAT A( J2 + 1 ) 0374 

. OLDSR=TEMPR 0375 

ULDSI=TEMfI 03/6 

J2=J2-IFP? 037/ 

IF( J2-J3)t75,675,670 03/8 

675 TEMPR=WP*SU,mf-0LDSR4JATA( J2I 0379 

TEMPI=V»I =i‘SUMI OjeO 

WORK{iO=TEMPF.-TEMPI 0381 

WORK IXCOM!) ) =TEMFP fTEMPI 036^ 

TEMPR = WR*5UMI-OLD5Ii-OATA( J2<-1) 0363 

TEMPI = WI «SU^1R 0384 

W0RK(K<-l)=TEVPR«-Tt;'1PI 0385 

WORKIKCOWJ + 1) =TEMPR-TEMPI 03 o6 

680 continue 0367 

IF(KN[N-1)685,685,686 0388 

685 WR=I«STPR<-1. 0389 

WI=WSTPI 0390 

GO TO 690 0391 
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686 TEM^'R = ^«(5 0392 

i^r=WP*WSTPR-K1=xwsTPI fWR 0393 

Wl = TEMPR*wSTPH-w!*WSrPR«-Wl • 0394 

690 TwOlNP = WPfWR 0393 

IF( ICASE-3)o92f 691,692 039o 

69L IF(IEPI-NP2I695,692,692 0397 

692 K=1 0398 

12MAX=I3*-NP2-NP1 0399 

DO 693 1 2=1 3, T2MAX,N.P1 0400 

DAT A( 12 I =i'!0PK( K I 0401 

DATA n2«- 1 )=WURK (Kfl) 0402 

693 K,= K+-2 0403 

GO T<^ 698 0404 

C 0405 

C complete a real transform in THt 1ST DIMENSION, N ODD, 8Y CG 0406 

C JUGATE symmetries AT EACH STAGE. 0407 

C ” 0408 

695 J3MAX=I 3+IFP2-NP1 0409 

00 697 J3= 13, J3MAX,NP1 0410 

J2MAX=J3 + I'P2- J2SIP 0411 

DO 697 J2= J3, J 2MAX,J2STP 0412 

J 1MAX=J2»J LRG2- IFP2 0413 

JlCNJ=J3fJ2MAX+J2STP-J2 0414 

DO 697 J1=J2, Jl!«AX,IFP2 0413 

K=1+J1-I3 0416 

OATA(Jl)=wnRK(K» 0417 

0ATA(JH-1) = WC’RK (K4-1) 0418 

IF( J1-J2 1697,697,096 0419 

696 DATA! JICNJ )=WGRK(K) 0420 

OATA( JICNJ+ l»=-Y.CPK(t9+l» 0421 

697 J1CNJ= J1CNJ-1FP2 ' 0422 

698 CONTINUE 0423 

IF=IFfl 0424 

IFP1=IFP2 0425 

IF( IFPl-NPi 1 7CC ,700,610 0426 

C 0427 

C COMPLETE A PEAL TRAN3FCRM IN THE 1ST UIMENSIJN, N EVEN, 8Y C 0428 

C JUGATE SY»-METRIES. 0429 

C ’ 0430 

70C GOTO <90C,8GC,9Cn,70l» ,ICASE 043i 

701 NHALF=N 0432 

N=N+N 0433 

THETA=-TWOPI/FLCAT (N ) 0434 

1F{ ISIGN >703,702,702 0433 

702 THFTA=-THETA 0436 

703 SINTH=SINITHETA/2.) 043/ 

WSTPR = -2.’<‘SINTH>>SINTH 0438 

WSTPI = S1 N( THETA ) 04j9 

WR = WSTPR<-1. 0440 

WI=WSTPI 0441 

IMIN=3 0442 
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JMIN=2*NHALF-1 
GO TO 725 
710 J=JMIN 

00 720 1=IMIN,NTCT,NP2 
SUMR=(0ATA( n<^DATAlJ> )/2o 
SUMI = (OATA( If 1 ) fCATA ( Jf lH/2. 

□ IFR= (DATA( I 1-0ATA( J ) 1 /2. 
OIFI=(OATA( If 1 )-OATA< Jf 1) >/2, 
TEMPR=WR*SIIMI f,V( I«OIFR 
TEMPI=WI*SUMI-WF«U1FR 
OATAC I I = SUMPf TENIPP. 

OATA( Ifl )=OIF I fTEMPt 
OATA( J)=SUMR-TEPPP 
OATA( Jfl )=-DIF IfTEMPI 
720 J=JfNP2 

IMIN=lMINf2 

JMIN=JMIN-2 

TEMPR=WP 

WP-WP*WSTPR-WI*V.STPI fWk 
Wl=TEMPP«'WSTP! f Vs I>s=wSTPP fV«l 
725 IF( IMIN- JMI N) 7 1C, /30,7^r0 

730 IF( ISIGN)731 ,74C, 740 

731 DO 735 I=IMIN,NT0T,N“2 
735 JATAdfl )=-DATA( Ifl) 

740 NP2=NP2fNP2 

N70T = MT0Tf fs!TOT 
J=NTOTf 1 
I MAX=NT0T/2f 1 
.745 IMlN=Ii'1AX-2«NHALF 

I=IMIN 
GO TO 755 

750 DATA(J)=DATAn » 

DATA! Jfl )=-L)ATA{ !♦!) 

755 I=If2 

J=J-2 

IF( I- IM4X ) 750, 7feC, 760 
76 0 DATA! J ) = OAT At I M IN )-UAT A ( I MINf 1 ) 
OATAIJfl )=C« 

IF! I-J > 770, 78C, 780 
765 DATA! J ) =DAT A( I ) 

DATA (Jfl )=DATA( Ifl) 

77G 1=1-2 

J = J-2 

IF( I-IMIN) 775, 775, 765 
775 DATA(J)=DATA(IMIN)fDATA(lMINfl) 
OATAIJfl )=0. 

IMAX= IMI N 
GO TO 745 

780 DATA! 1 )=04TA( 1 ) fDATAi 2) 

0ATA(2)=0. 

GO TO 900 


0443 

0444 
04-^3 
0440 
0447 
04'»6 

0449 

0450 

0451 
sJ45Z 

0453 

0454 
j455 
0456 
045 7 

0458 

0459 
046U 
0461 

0402 
0463 
04t>4 

0403 
0 466 
04o 7 

0468 

0469 
04 70 
0471 
047z 
04 73 
04 7 4 
04/5 
04 7 6 

0477 

0478 
04 79 
O48o 
0461 
048 2 

048 j 
-,.04b‘t 

o465 
C 466 
04 o 7 
0488 
0469 
0490 

049 1 
0492 
049 3 
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C , 

C CCMPUtTE A REAL TRANSFORM FOR THE 2NU OR 3RO UlMENSION BY 04V5 

C CONJUGATE SYMMETRIES. 04V6 

C .0457 

800 IF( IIRNG-NPU 8C5,900,900 049e 

805 00 860 I 3=1 ,NTCT,NP2 0499 

I2MAX=I3«-NP2-NPl 05CC 

00 860 I 2=1 3, I2MAX.NPI 0501 

IMIN«I2fIlPNG 050*; 

IMAX=I2fNPl-2 05u3 

JMAX = 2*I 3*-NPl-IMlM 05C4 

IF( 12-13 )320,820,810 0505 

810 JMAX = JMAX«-NP2 05C6 

820 IF( IOIM-2) 850, £50,830 0507 

830 J = JMAX<-NPO C5Lio 

00 840 I = 1MIN, IMAX,2 05u9 

OATA( I ) = OATA( J ) 051C 

OATA( I4-1)=-DAT A( J + 1) 0511 

840 J=J-2 0512 

850 J=JMAX C513 

00 860 I = IMIN, IKAX.NPO 0514 

OATA( M = OATA{ J ) 0515 

DATAIUl )=-OATA(J + l) 05io 

860 . J=J-NPO 0517 

C 0518 

C END OF LOOP ON EACH UlMENSION 0519 

C • 0520 

900 NPO=NPl 0521 

NP1=NP2 0522 

910 NPREV=N 0523 

920 RETURN 0524 

END 0525 

OVERLAYIPATS, 1 ,C) 

PROGRAM P.EAOIM 

COMMON/BLK l/STAPTT.l TFMT.NBLK, I P0w2 , NCH , NPR I N T , I PLOT A, IPLUTC.OFFSC 
1AL,0ELTAT, SN,NR SKIP, LAP.NC ROSS, I CROSS (2,20) , NCHP , YLA BE L I 2 I , IWINUOW 
1,P1 ,F2, ITYPESP ,hCCN,NFSKIP,IFF,LAGl,LAG2 
C0MM0N/3LK2/ICHI 14) ,CHSL'M( 14) ,N0FFI14) ,CHSUMSO( 14) , SIGMA ( 14) ,R.MS( 1 
14) ,MEANl 14 ) ,SCALFAC( 14) , ChSUM 1 U 4 ) , T R AC K ( It ) , 1 CHAN I 14),1FILTER114> 
COMMON/BLK 3/ NPT ,TMAA ,NPT02,NSPCT ,OELF ,N6 4, NPT Cl 28 
1,IN2ER0,NREA0,NPT0T 
C0MM0N/3LK4/NMAX,NCHMAX ,N8CMAX 

C0MM0N/BLK5/PCTC,NBINS,0MAX( 14 ) , OMIN (14 ) , OB I N ( 14 ) , B IN S ( 1 GO , 14 ) 

IfCHISOC 

C0MM0N/BLK8/I AUTQSPI 14 ) , I AUTOCOI 14) , ICRS P ( 20 ) , I CRC(JR(2C) , I TRA (20) , 
1IC0H(20) 

COMMON/BLK 9 /NF I L TP , F H EQ F ( 50 ) , WGHT F ( 50 ) 

NAMEL 1ST /I NPUT/ ITFMT , NF SK I P , NRSK I P, SN,UELTAT , ST ART T , OF F SC A L , NCH , 
ISCALFAC, NPTOT, NPEAO, lAUTOSP , I AlJT OCO, NCR OS S , I CROSS, I CR5 P ) IC RC Oh , 
2ITRA,IC0H,LAP, I hlMDOW, I T Y P E S P , NPR I NT , I PLUTA, 1 PL OTC , F 1 , F 2 , L AG 1 , L AG2 
3 ,PCTC,NBINS,DMAX,UM1N, IN2tR0 ,NFILTP, FREUF,wGHTF,IF ILT tK 
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DATA NFILTP,HPECF,WGHTF, IHILTcR/0,100^0. 

DATA IAU’’OSP, I A L TOCO, ICRSP, ICKCUK, ITkA, ICOH/106*0/ 

DATA IFF,LAG1,LAG^/1,0,0/ 

DATA NFS KIP/0/ 

DATA ITFyiT, STARTT,ICH,NRSKIP, SC A L F AC , OF F SC At , L A P , NCROS S , I CROSS 
INPR IMT, I PLOT A, IPLOTC, I W IN DOW , Fi , F 2, I T YP E SP/ 

22 ,0* , IS*C, ,l»E*'t,42*0, lOu, »-<• ,2 00 00* , 2 / 

DATA TD,CHAN/2HI0,4HrHAM/ 

DATA PCTC/90./ 

DATA lNZEPn,NBINS,0yA.X,D,KlN/2*O,26*O./ 

READ AND PRINT NAMELIST INPUT 

RfcAO INPUT 
IF(£0F,5) 1,2 

1 STOP 

READ AND PRINT FOPMATTFu CARO INPUT 

2 PRINT INPUT 
READ 5002, YLA3EL 
READ 5002, (TRA CK ( I ) , 1= I, NCH» 

50C2 FOpyATIBAlC ( 

PRINT 5001, YL AbEL 

PRINT 5003, ( I , TRACM I ) , I =1 ,NCH> 

5001 FOP.yATI//^ CASE ID*5X 1H*2A 10, iw* ) 

5C03 FORMAT!//* CHANNEL I D*/ 3X *NU, ♦/ ( I 5 , 4 X A1 0 ) ) 

PRINT 500A 
500A FORMAT!*!*) 

r. 

C CHECK INPUT DATA, PRINT ERROR MESSAGES 

C 

NPT=NRE AD 
nrlk=nptot/npt 
IF!NCH,LE.NCHMAX) go Tq Id 
PRIN'^ 103, NCHMAX 
STOP 103 

103 FORMAT(///* GT*I3,*, PROGRAM. wILL NUT READ TAPE CORRECTLY. 

1 terminated*) 

101 CONTINUE 

IF!NFSKIP) 310,210,311 
ill DO 312 I FSK IP= 1 ,NFSk IP 
313 READ!1) SKIPRFC 
IF!EOF,l) 312,313 
312 CONTINUE i PRINT 90A , NFSKIP 

904 F0RMAT!//6H * * *,!5,20H FILES SKIPPED * * *) 

310 IF! ITYPESP'.LT.31 GO TO 9C 
NCR0SS=0 
DO 91 1= 1,NCH 

IF!IAUT0CD! II .GT.C) lAUTOCO! I )=0 
91 CONTINUE 


J Ob 
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PRIMT 92 

92 "FORMAT ( //=)■ AMPLITUDE SPECTRUM OPTION CH0SEN'>/=!‘ ONLY AUTO SPECTRUM 
IW ILL BE C ALCUl ATED*) 

90 CONTINUE 
IC0RR=0 

DO L16 1=1, NCH 

1F( lAUTOCOI I) .EC.O) GO TO U6 
IC0RR=1 

IF( lAUTQSPI n.EC.O) lAUTOSPII »=-l 

116 CONTINUE 

IF(NCROSS) 118,116,108. 

108 UG 110 I=1,NCR0SS 

J1=ICP0SS( 1,1 I $ J2= ICROSSI 2, n 
IF( I AUTOSP (J1 ) . EQ.O) I AUTOS P« Jl)=-1 
TF(IAUTOSP(J21.EO.O) 1AUT0SP(J2)=-1 
IF( ICRCDR ( I ).NE.G» GO TO 117 
1F( ITRAI 1 ) .NE.O » GO TO 21 1 
IF( ICOH( I I.NE.OI GO TO 211 

GO TO no 

117 IC0RP=1 

211 IF( ICRSP ( I ) .EQ.C » ICRSP(II=-1 
no CONTINUE 

118 1F( ICORft.EO.l.AND.IN/ERQ.EO.O) PRINT 90 7 

907 FORMAT(//ir 7H YOU MAY HAVE CIRCULAR ERROR IN YUUR CORRELATIO 

INS BECAUSE YOU HAVE NUT ASKED FUR ZERO INSERTION * * *1 
IF( IMZERO.EQ.O) GO TO 109 
NPT = 2<'NPT 
PRINT 901, NPT 

901 FORMAT(//'C ZERO INSERTION INCLUDEU FOR ALL CHANNELS. BLOCK SIZE IS 

nm 

■ 109 do 3f> IK = 1 ,NCH 
N0FF(IKI=0 

CHSUMI IK )=CH5UMSQUK )=0. 

CHSUMl ( I Kl =0. 

35 CONTINUE 

IPOW2=0 $ NTEMP=NPT 

203 IF(NTEMP-l) 200,200,202 

202 IP0W2=IP0W2*-1 $ NTEMP = NTEMP/2 i GO TO 203 

200 IF< 2*<‘IP0W2.NE.NPT) GO TO 201 
PRINT 207 

207 FORMAT!//* BLOCK SIZE IS A POWER OF TwO. FAST FOURIER TRANSFORM WI 
ILL 3E USED*) 

IF( NPT.LE.NMAX ) GO TO 205 
PRINT 206 $ STOP 207 

201 PRINT 204 

204 FORMAT!//* BLOCK SIZE NOT A POWER OF TWO. SLOW FOURIER TRANSFORM W 
lILL BE USED*) 

IF !NPT, LE.NMAX ) GO TO 2C5 
■PRINT 206 

206 FORMAT!//* BLOCK SIZE TOO LARGE FOR DIMENSIONS PROVIDED. JOB TcRMI 
INATED*) 


86 



o o o 


APPENDIX H 


STOP 2Co 
2C5 CONTI. NiUE 

NPT02=NPT/2. 

NSPCT = NPLm =NPTC2 i I F ( NPP INT .GT .NPT02 > NPR1NT=NPT02 

TMAX = NPT#nELT .'-T 

OELF=U/TMAX 

NPT0128=NPT/128 

IF( INZERG.EU.O.GR.LAP.EO.C ) GO TO ill 

PRINT 902 

LAP=0 

9J2 FORMAT{//>:= N.' .SO PERCENT OVERLAP ON ZERO UmSEkTION RuNS«/* INPUT 
LOATA ALTERED, LAP=0 
111 CONTINUE 

IF(LAP.NE.O) NbLK = 2«Ne-LK - 1 
NCHP=0 

DO 80 1=1, NCH 

IF( IAIIT0SP( I) I 8l,UG,8l 

81 NCHP=NCHPvl 
ICHANfNCHP 1=1 

80 CONTINUE 

IF(NCHP) 82,82,83 

82 PRINT. 84 

84 FORMAT!//-!-' INPUT INDICATES NO CHANNELS TO oE PROCESSED, CASE ENDED 
1*) 

STOP 101 

83 CONTINUE 

IF(NCHP«NlU.K.LE.N8CKaX) GO TO 3C6 
PRINT 30 7, NCHP ,N8LK ,NBCMAX 

307 FORMAT!//'!' NU . OF CHANNELS TO BE PROCESSED !NCHP» =<‘IS/>i' NO, OF BL 
lOCKS (NRLKI =*I5//33H NCHP * NBLK GREATER THAN NDCMAX= , I!>/ / * EXECU 
2TI0N ENDED. CHANGE DIMENSIONS TO FIT YOUR CASE AND RERUN-*-) 

STOP C7 
306 CONTINUE 

N64=NCHP«64 

C 

C COMPUTE ACCtJRACY MEASUREMENT OF SPECTRAL ESTIMATORS 

C ' 

IFILAp.EQ.O) NDCF=2«NHLK i If ILAP.NE.U) NDuF = 1 . 6364^ ! NBL K- 1 ) 

CALL CSQ ! ! lOG . + PCT C) / 2 . , NDOF , BU, B L, I C ODE ) 

BL=1C0./BL t BU = 1C0. /CU 
PRINT 903, PCTC,BI.,HU 

903 FORMAT ! ///58H * * ACCURACY MEASUREMENT OF SPECTRAL ESTIMATORS 

i* * -S.//-X .ASSUMING normality OF DATA, USER. CAN BE<-F6.0,* CERTAIN*/ 
2/* THAT THE SPECTRA! ESTIMATE IS wlThlN THE BOUNDS OF -»//F5.0* PE 
3RCENT ANO-i'FB.O,* PERCENT OF THF TRUE SMOOTHED SPECTRUM*) 

COMPUTE CRITICAL VALUE OF CHISQUAK6 FOR NORMALITY TEST 

IFINBINS.GT.ICO ) NRI,NS=100 
IF!NRINS.FU.O) GO TO 113 
DO 112 1=1 ,MCHP 
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<=ICHAN( I ) 

bBI:Ml 11= (DM AX( K )-i)MIN(K > ) /NBINS 
DO 112 J=1,NBTNS 

112 fiINS( J, I )=0. 

NBlNSM3 = N'UNS-3 

CALL CSO(PCTC,NBINSM3,BL,BU,ICUOfc ) 

CHIS0C=NBINSM3/BL 

113 CONTINUE 

IF( NEILTPolO.O I GO 10 114 
IE(Nf (LTP.LF.BCI GO TO 115 
PRINT 905 

905 FORMAT ( ///55H * ★ INPUT ERROR, NFILTP GT 50. EXECUTION ENDED ♦ 
1 

STOP 06 
115 PRINT 90b 

906 FORMAK //126H * * ^ SPECTRAL FILTERING OPTION SELECTED. CORRELATIU 
INS WILL BE CALCULATED FROM FILTERED SPECTRA FOR CHANNELS SPECIFIED 
2 . * 

114 CONTINUE 

COMPUTE DATA WINDOW CORRECTION FACTUR FOR SPECTRA 

CAPT=ORE AU#UEL TAT $ I W I NOP 1 = I W I NUUW + 1 
GO TO ( 300,301,302,303) .IHINDPI 

300 WCQN=CAPT $ GO TO 304 

301 WC0N = CAPT>!<. 375 $ GO TO 3C4 

302 WCON=CAPT*. 7136727029 t GO TO 304 

303 WCQN=CAPT*. 2696428571 

304 CONTINUE 

CORRECT PLOT LIMITS 

IF(FULT.O.) Fl=Co 
FMAX=NPT02*OELF 
IF(F2.GT»FMAX) F2=FMAX 
IFILAGI.LT.-NSPCTI LAGl=-NSPCT 
, IF(LAG2.GF.NSPCT) L«G2 = NSPCT-1 
IF(LAC2.GE.LAG1 ) GO TO 305 
L1=LAG2 t LAG2=LAG1 i LAG1=LL 

305 CONTINUE 
RETURN 
END 

SUBROUTINE CS 0 ( P , N , B L , BU , IC006) 

EXTERNAL FUNC 
C0MM0N/'’RD3F/ A,t 
A=N/2. 

0=P/100. 

C FIND UPPER BOUND 

CALL I TR 2 ( CHI S g ,.CQ1 , iOOC. , lo ,FUNC, I.e-6,l„E-6, ICO, I CODE ) 
IE(ICnuE-3) 1,1,2 
2 BU=1.EH2 
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GO TO 3 
1 BU=N/CHISO 

3 y=l.-Q 

FIND LOWE« BOUNC 

CALL 1 TP.2(CHI SOt.COl , 1000. ,1. ,FUNC, i.E-6 , 1. E -6 , ICG , I CODE ) 
IF(IC0DE-3) 4,4,5 

5 BL=0. 

GO TO 6 

4 8L=N/CHI SO 

6 RETURN 
END 

FUNCTION FUNCICH/SQ) 

COMMON/PROBF/A ,C 
DIMENSION H60( 15) 

DATA H60/'”.01lB,"“.0o67,~«0033 ,~.‘j010,.0CCi,2’^.u0CB,.G0w2,“'.0C0^,, 
1-.0006, - .00 G5 ,. C0G2, .0017, .0043, . 0062/ 

IF( A-15) 1,1,2 

1 X=CHlS0/2. 

FUNC = 0-GAMMF( A,XI/GAMMf (A,0. ) 

GO TO 3 

2 X=( (CHISO/ (2.<=A ) )*>:=( 1./3. )-( i.-1.7(S,.*A ) ) l/SQRT( l./(0.*A I ) 

IX=X/.5 $ !F(IX«.5.GT.X) IX=IX-li IX=IX+8 i IF(IX) o,6,/ 

7 I F ( I X- 1 5 I a , 8 , 6 

8 IF( IX.LT. 1 ) GO TO 6 
XREM = X-( IX-b)«. 5 
IF(XP.EM) 4,9,10 

9 HNU=(30./A)*H60(IX) 

GO TO 11 

10 HNU=(30.7A)*(H6C(IX)+XREM*(HC>0( IXfl)-HoO(IX) l/o5) 

GO TO 11 

6 HNU=0. 

11 CONTINUE 
X=X<-HNU 
AX=ABS(X) 

X2=AX>i‘AX $ X3 = AX>i‘X2 $ X4=AX*X3 

PMINUS1 = -. 5*( 1. ♦. 196 854’^AXt. 1 1 5 1 94*X2<-. 000344X=X 3 *■. 01 9527'»'X4 / ♦♦( -4. 
1) 

1F{ X) 4, 5,5 

4 FUNC.=Q-( 1.<-PM!NLS1) 

GO TO 3 

5 FUNC = y<-PMlNUSl 

3 RETURN 
END 

0VERLAY(PATS,2,G) 

PROGRAM BLOCKS 
COMMON CMAIN( 1 » 

CCMM0M/BLK1/STARTT,I TFMT,NbLK, I P 0W2 , NCH , NPR 1 NT , I PLOT A, IPlOTc,UFFSC 
lAL, DELTA T, SN , NR SK I P , L AP , NC ROS S , IC ROSS (2 ,20) , NCHP , Y LABEL < 2 ) , IWInUOW 
1,F1 ,F2, I TYPESP, V>CCN, NFSKIP , IFF,lAG1, LAG2 
COMMON/ BLK2/ 1 CHI 14) ,CHSUM( 14) ,NJFF < 14) , CH SUM SOI 14 ) , S I GM A ( 1 4 ) , RMS ( 1 
14),MEAN( 14) ,SCALFAC( 14) ,CHSUM1(14),TKACK(14) , ICHANI 14) 
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C,0MMQN/BLK3/NPT , T MAX , NP T 02 i NSPCT , UELf- ,N6 4,NPT0i2b 
1 , INZ6R0, NREAD 

COMMON/ BLKD/pr.TC,N3 IMS , DM A X( 1 ^ , DM I N ( 1 ^ , DB I N ( 1 4 ) , B 1 NS ( 1 cO , 1 4 M 

l.CHISUC 

C0MM0N/BLK6/I S AVE64, IR I , IXPLOT, I DATA , I Z , I SPcCT 
C0MMn.N/BLK7/I T ( 11) 

C0MM0N/BLK8/I AUTQSPM 14 » , I AUfOCGl 1 4 ) , I CR S P ( 2 0 ) , ICRCORl 20 » ,1 TKA( 20) . 
HC0H(20) 

00 1 NB=1,NBLK 

CALL P.CADTPE (NB tCMAIf'U I DATA ) ,CMA INI ISAVEo4) ) 

CALL TRAN(NB,CMAlfM( I Zl , CMAIN( ISAVE04) ,CMA1N( I SPECn ) 

1 CONTINUE 
RETURN 
END 

SU8R0UT1 NE RE ADTPE (NB, DATA, SAVEo4 ) 

DIMENSION OATA( 1),SAVE64( n 

COMMON/BLKl /STARTT,! TFMT.NBLK, ITQw2 tNCH , IMPRINT, I PLUT A , I PLOT C , OFF SC 
1AL,DELTAT,SN,NREKIP, LAP.NCROSS, ICROSS (2 , 20 ) , NCHP , ¥ L A BE L ( 2 ) , IWINDOtM 
l.Fl,F2, I TYPESP,RC0N,NFSKIP,IFF,LAG1,LAG2 
COMMON/ BLK2/ 1 CH (14 ) , C HS UM ( 14 ) ,N OFFC 14 ) , CHSUhSO( 14 ) , SIGmA( 14) ,RMS( 1 
14) , MEAN( 14) , Sr ALFACI 14 1 , C HSUM 1 ( 1 4 ) , T R A C K ( 1 4 ) , 1 CH AN ( 1 4 ) 
C0MM0N/BLK3/NPT ,TMAX ,NPT02,NSPCT ,OELF ,Nb4,NPTC12B 
1 , INZERO, NP'E AO 

C0MM0N/BLK7/NN , I WO ,K CH , NFR, 1 1 , KK , N1 , J J , I REC, NR E C , 1 L OC 
IF( I TFMT .E0.3 ) GO TO 4000 
IF(NB-l) 12,12,11 
12 GO TO ( 3001,3010, ITFMT 
C AOTRAN FORMAT PARAMETERS SET UP 

3001 IF(NPSKIP) 105,105,106 

106 00 107 I =1 ,NRSK IP 

107 REAO(l) 

105 CONTINUE 

IF(NCH.GT. 10) GC TQ3002 

KCH=20 

NFR=25 

GO TO 3003 

3002 1F(NCH.GT.20) GC TO 3004 
KCH=30 

NFR= 17 
GO TO 3003 

3004 IF(NCH,GT,30) GO TO 3005 
KCH=40 

NFR=12 
GO TO 3003 

3005 IF( NCH,GT.4C) GC TO 3006 
KCH=50 

NFR=10 
GO TO 3003 

3006 1F( NCH.GT. IOC ) GO TO 3G0Z 
KCH=110 

NFR = 4 
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GO TO 3003 

3007 PRINT 300H 

3008 FORMATC//# NCH GT 100 NOT ALLOwfcO<=) 

STOP 02 

3003 IW0=9 

NN=KCH«NFP 
GO TO 3011 

C AOTRAN interface FORMAT 10 RECORD 

3010 REAOm KFV,NN,IwO,KCri,NFR, lO'i, I02iTSN 
PEAOm 

READ! 1) 

1F( TSN.EO. SN) GO ■^0 301 1 
PRINT 3012 

3012 TAPE Ni^T PPSITIONEO AT ID RECORD FOR DESIkEd SN*) 
STOP 3012 

3011 IF(NRSKIP) 3009,3009,108 
lOB on 109 I=1,NKSKIP 

1C9 READ! 1 ) 

FIND STARriNO TIME ON TAPE 

3009 CONTINUE 

Gn TO (3013,3014), ITFMT 

3013 REAO(l) (DATA ( I ), 1=1 ,NN I 
I F ( t OF , 1 ) 3 , 2 

2 IF ( SN. Eg.OATA ( 2 I ) .GO TO 3C15 
PRINT 3016 

3016 FORMAT!//* TAPE NCT POSITIONED AT UtSlREU SN*) 

STOP 04 

3014 READ(l) KFY,NM, (DATA! I > , I = 1,NN) 

I'F(E0F,1) 3,3C15 

3015 DO 1 J = 1 ,NFR 

I 1= ( J-1 ) *KCHf 1 WD»-1 
JJ=J 

IF(DATA( I' )-STARTT) 1,4,4 
1 CONTINUE 
GO TO 3009 

3 PRINT 9C0, STARTT 

900 FORMAT!//* STARTING T I ME *(i 1 2 , 5 , * NOT FOUND ON TAPE*) 

STOP 

4 KK=0 
N1 = J J 

START READING DATA TO BE PROCESSED 
GO TO 13 

il N1=JJ+1 $ KK=C 

IF(NI.GT.MFR) go to 15 
BACKSPACE 1 
GO TO 16 

15 Nl=l t I I=l WUF 1 
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16 IF( ITt^MT.Pg.l » FE/10(1) (UaTA(n,I=l»NN>- 

IF( I THMT .fcO.2 ) RFiO( L) KE Y , NN , ( li A T A ( I ) , I = 1 , NN ) 

IF(E0F,1I 30o3f3066 
3066 CONTIi^UE 

IF(LAP) 13,13,3062 
206 ILUC=I1-1 

IREC=IREC1-1 
GO TO 14 

13 IREC=0 

14 ILOC=C 
10 CONTINUE 

DO 5 J=N1,NFR 
KK=KK+1 

IF(KK-NREAD) 6,6,7 • 

6 IL0C=IL0C«-1 

IJ= ( ILOC-n *NCHP 
DO 30 IKK=l,NCHP 
IK=ICHAN(TKK) 

X = OATA( I UIKJ =»SCALFAC( IK I 
IF(ABS( Xl-OFFSCAL) 33,34,34 

33 CHSUM(IK»=CHSUM(IK|+X 
SAVE64( I JMKK ) = X 

GO TO 30 

34 SAVE64( IJfIKK » = CHSUM ( IK »/ (KK-l4-(NB-l)*NRtAD» 

CHSUMI IK » =CHSUM( IK)* ( KK ♦■ ( NB-1 ) *NR EAD > / ( < KK* ( NB- 1 ) *NRE AD ) +1. ) 

NOFF ( IK ) =NOFE ( IK) n 

30 CONTINUE 
IF(IL0C-64) 31,32,32 

32 ILQC=0 $ IREC=IREC*1 

CALL WRITMS(9,SAVE64,N64, IREC) 

31 CONTINUE 
J J=J 

IT=II*KCH 
5 CONTINUE 

IF(KK-NPT) 8,7,7 
8 CONTINUE 

IF( ITFMT .EO. 1 ) REAU(l) ( U A T A ( I ) , I =1 , NN ) 

IF( I7FMT.E0.2 » REAO(l) KE Y , NN , (DATA ( 1 ) , 1 = 1 , NN ) 

IF(E0F,1) 3063,3064 

3063 PRINT 3065 
STOP 3065 

3065, FORMAT!//* END GF FILE ENCOUNTERED BEFORE. NPTOT POINTS READ, EXECU 
ITION ENOEO*) 

3064 CONTINUE 
Nl=l 

1I=IWD+1 
GO TO 10 

ONE BLOCK REAL, READY TO BE PROCESSED BY FFT 
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IF( RnC.EQ.O) RETUKN 
IREC=IREC«-l 

CALL WRI TMS(9,SAVEb4,lLQC*NCHP,IREC) 

RETURN 

RECIN FOR'^AT 

4000 .NREC = NRE AO/64 $ NCHP2 = NCHf2 $ Il=i 
1=0 

IF(NH-1» 4GC1 ,4CCl,3029 

4001 IF(NRSKIP.EO.O) GO T0'30t)0 
on 3051 I RSKI P= 1 ,N°SK IP 

3051 CALL REC IN ( 1, 2 ,NIN ,QATA , 1 ,NCHP^, 1 ) 

PRINT 3057, NPSKiP 

3052 FORMAT(//I 10,« RECUPOS SKiPPEDi--) 

3050 CONTIMUE 

1 = 1^1 

CALL REC IN ( 1 ,2 f KN.UATA , 1 ,.\CHP2, 1 ) 

IF(£GF,1» 3,3L27 

3027 IFIDATAI 21-STARTT) 3 0 50 , J C2 2 , 302 2 

3022 IF(i)ATA( D.EQ.SN) GO TO 3L35 
PRINT 3016 

STOP 05 
3035 !1=2 

PRINT 901, I 

901 FORMAT (//-K STARTT FuUND AT RECCRl)«=I5) 

KK=1 

on 3023 , 1KK = 1 ,NCHP 
IJ=ICHAN(IKK) 

XX=OATA( IJ«-2) «SCALFAC( I J ) 

IF( AflSI XX) -OFFSCAL ) 30 lO , 30 31 , 30 3 1 
3031 XX=0. $ NUFFI I J )=NOFF( I JH-1 
3030 CHSUMI 14 ) = CHSIIM( IJ )fXX 

3023 SAVE64( 1 KK ) =XX 
IREC1=1 

NREC = NRE AO/64 t GC:' ;3 3C6C 
3029 IF(LAP) 3061,3061,3062 

3061 KK = 0 t 11 = 1 t I REC 1=1 % NREC=NRE AO/o-t $ GO TO 3060 

3062 KK=0 
IWRITE=C 
IPEC=NPTa2/64 
NREC=NRE Ai)/64 
IL0C=(NPTn2-IRE C*66) «NIC HP 

IREC=IREC+1 $ NW0R0S=N64 t NL EFT =No4- I L JC 
4018 IF ( IRFC-NREC) 4 C 16 , 40 lo , ^0 1 7 
4017 NWORDS= ( NRE AD-NREC=!=t.4 ) -^NC HP 
NLEFT=MWORL)S-ILOC 

4016 CALL READMS(9,SAVE64 ,NWORDS, IREC ) 

IREC= IREC^l 

KK=KKfNLEFT/NCHP 

IF( ILOC. EQ.O) GO TO ^014 
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OG 4006 I=l,NLtfT 

4006 SAVE64( I ) = SAVE64( IHLOC ) 

IF(KK-NPT02) 4007,4012,4013 

4007 IF ( IREC-NREC) 4004,4009,4010 
4C10 NWOROS= ( NRE AD-NREC*64 )>!>NCHP 

IF( NWORDS. LT.. ILCC) ILOC=NWORDS 
4009 CALL P.EAOMS (9 , SAVE64( Nl EFT + 1 I , ILOC, 1 PEG I 
KK=KK«-ILnC/NCHP 

4014 IF(KK-NPTn2» 4011,4020,4021 

4011 IWR ITE= IWR I TE «■ I 

CALL WRI TIN(9 , SAVE64 ,N64, I WRI TE ) 

GO TO 40 13 

4012 IF( NLEFT.lt. N64) GO TO 4019 
4022 IWRITE= IWR ITE*-1 

CALL WRI TIN(9 , SAVE64,N64, IWR ITE) 

11=1 

4015 IREC 1= IWRI TE+ 1 
KK=NPT02 

GU TO 3060 

4013 I1 = NLEFT/NCHP-(KK-NPT02) «•! 

GO TO 4015 

4019 Il=NLEFT/NCHPf 1 
GO TO 401-5 

4020 IFt IL0CfNLEFT.EC.N64) GO TO 4022 
Il=( ILOCfNLEFT ) /NCHf f 1 

GO TO 4015 

4021 Il=(NLEFTf ILOC ) /NCHP - ( KK-NPT02 ) f 1 
GO TO 40 15 

3060 IF( ITFMT.lt. 3) GOTO 206 

IF( IKECl.GT.NREC) GO TO 3046 
DO 3024 IREC= IREC1,MREC 
00 3025 IL0C=I1,64 
CALL REC INI 1, 2 ,NN,OATA, 1 ,NCHP2, 1 ) 
IF(E0F,1) 3063,3034 
3034 KK=KKfl 

IJ=( ILOC-1 )*NCHP 
00 3026 IKK=1,NCFP 
IK=ICHAN(IKK) 

XX=DATA( IKf2) «SCALFAC ( IK ) 
IF(ABS{XX)-OFFSCAL) 303 2, 303 3 , 3C 3 5 
3 033 XX=CHSUM ( I K » / ( KK-i f I NB- I ) UNREAL)) 

NOFF ( IK ) =NGFF ( 1 K) t 1 
3032 CHSUMI IK ) = CHSGM{ IK)fXX 
3026 SAVE64I I Jf IKK. ) =XX 
3025 CONTINUE 
11=1 

3024 CALL WRI TMSI9 , S AVE64 ,N64, IREC ) 
IFINREC*64.E0.NREAU) RETURN 
3046 NLEFT=NREAD-64*NREC 
NREC=NRECfl 
00 3042 ILGC= I 1 ,NLEFT 


94 



r» o r> 


APPENDIX H 


CALL PEC INI li 2 ,NN,I)ATA, 1 ,AjCHP2, 1 » 

IFtEOE,!) 3063,30^1 

3041 KK=KKfl 
IJ=( UDC-1 

no 3042 lKK=l,NChE 
lK=ICriAN( IKK) 

XX=DATA( IK *-2) • 

If ( XX-OFFSC AL I 3043,3044,3044 
3044 XX=CHSUM( IK)/(KK-1<-(.‘J8-1)*MREA0) 

N0FF< IK) =IX1FF ( IK)+1 
3043 CHSUM( IK )=CHSUM ( IK )<-XX 

3042 SAVE64( I J+ IKK ) =XX 
CALL wPI TMS (4 , S AVF64 ,NLEFT*KiCHF,i\REC ) 

RFTUPf 
END 

SUBROUT I iMF TPAN(NB,Z ,SAVL64,SPECT) 

OIMENSION Z ( 1 ) , SAVE64( 1 ) , SP£CT( 1 ) 

COMPLEX Z 

COMMflN/aLKl /START! ,ITFMT,NBLK, I POw2 , KiCH , (X PR I (M T , I PLOT A, IPLUTC,OFFSC 
IAL,f)FLTA T , SN,!'JR SK I P, L AP , N C ROS S , I C PUS S ( 2 , 2 C ) , UCH P , Y L AO E L ( 2 ) , IW INIJOK 
L,F1,F2,ITYPESP 

C0MMQN/BLK2/K H ( 14) , CHSUM( 14) , MLF F ( 1h ) , C HSUM S 0 { 14) , SIGMA 1 14) ,RMS( 1 
14) ,MFAM( 14) ,SCALFAC( 14) ,CriSUMl(14),TRACKI14) , iCHANl 14) 

COMMON/RLK 3/NPT ,TMAX ,!^P TO 2 , tMSPCT , DEL F , N64 , NP TO'l 20 
1 , INZERO.NREAU 

CC!MMnN/BLK5/PCTC,KO.INS,OMAX(14),OMIM 14) , Ob 1 G { 14) ,BIUS( 100,14) 

00 19 JJ-1,NChP 

READ INPUT DATA RCOCK FROM RANOGM ACCESS FILE 

J=ICHAN(JJ) 

NREC=0 
IK=0 

IF(NRbAi).LT,o4 ) CO TO 21 
00 20 I = 64,iNREAC,64 
NREC=^NREC<-1 

CALL REA0MS(9,SAVE64,N64,NkEC ) 

DO 20 IJ=1,64 
IK=I-64KIJ 
IL= I lJ-1 )*NCHP ♦•J J 

20 Z( lK) = SAVto4( ID 
!F(NP.EC«64oEQ<,NFEAn) GQ TO 10 

21 MLEFT=NREAO-Ik 
NREC=NREC+ l 

CALL REA0MS(9, SA\,£64,NLEFT^«NCMP,rgRtC) 

DO 11 !J=l,NLKFT 
IL=( lJ-1 )*MCHPkjJ 
11 Z( IJK1K)=SAVE04 (IL ) 

10 CONTIUUE 

c 

C COMPU’^E BLOCK MEAN AND COUNTS FUR HISTOGRAM 
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0LKMF:AN=O. 

DO It I=l,MRfcAO 
CHSUMSO( J ) =CHSUMSQ( J ) fZ ( I 
2.2 BLKMEAN=BLKMEAN + Z( n 

CHSUMK J ) = RLKN|EAN*-CHSUMU 
BLKMEAN=BLKME AN/NkhAO 
IE ( NBINS .EO.O ) GO TO 40 
DMN=DMIN(J) $ DEL=UaiN(JJ> 

DU 41 I=1,NBEA0 
IRIN=(Z( I l-OMNI /DEL+l 
IE( IBIN.LT.l) IBIM=1 
IF( IBIN.GT, NBINS) IBIN=NbINS 
41 aiNS(IBIN,JJ)=BINS(IBlN,JJ)»l 
40 CONTINUE 

DO i!3 I = unread 

23 Z( 1 1=Z( I )-BLKMEAN 

APPLY DATA i^IMDCW 

IF( IWINOQW.EQ.D 1 GO TO 33 
GO TO (30,31,32), IWINDOW 

30 CALL HANNING! Z ,NREAD) 

GO TO 33 

31 CALL HAMMING(Z ,NPEAD) 

GO TO 33 

32 CALL PARZEN(Z,NREAO) 

33 CONTINUE 

insert zeros 

IF( INZERO. EO.O ) GO TO 2s 
DO 3t) I = l,NPT02 
35 Z( I <-NPT02) =C. 

COMPUTE FOURIER TRANSEOPM AND STORE ON RANDOM ACCESS FILE 

24 CALL EOURT! Z,NPT, 1,1,0, SPECT) 

IJ=(N3-1 )*NCHP«-JJ 

CALL WRI TMb (ti , Z ,NFT, 1 J I 
19 CONTINUE 
RETURN 
END 

SUBROUTINE HANN ING( Z , NPT ) 

COMPLEX Z( 1 ) 

DATA PI / 3. 141 5926535B90/ 

DO 1 1=1, NPT 

0 = .5>^'(1. -COS( 2. I-l. ) / NPT ) ) 

1 Z( I )=Z( I )-<D 
RETURN 
END 
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SUBROUTINE HAE'M I NG( Z , NPT ) 

COMPLEX Z(l) * CATA P I / 3, 1 4 139265 3S896/ 

PION = PI/NPT $ NF2=NPT*-2 
DO 1 1=1, NPT 
T={ I M-NP2) «^PTON 
0=.54f.46*CnS(T) 

1 zn ) = Z( T )=>!) 

RETURN 

END 

SUBROUT I NE PARZEMZ.NPT) 

COMPLEX Z( I ) 

TM=NPT/2 % TMM1=tm-1 
TM02=TM/2. 

DO 1 1=1, NPT 

j=1-TMM1 

IJ=IABS(J» 

1E( IJ-TM02) 2,2,3 

2 U=l.-6.^= J* J*( l.-l J> 

GO TO I 

3 D=2.4( 1. -I J |«^<3 

1 Z{ I ) = Z( I 
RETURN 
END 

OVERLAY! PATS, 3 ,C ) 

PROGRAM AUTOSP 
COMMON CMAIN(l) 

COMMQN/BLK l/STARTr,I TFMT,N6LK, I P0«2 , NC H , NPR I NT , I PLOT A , 1 PLUTO, OF ESC 
1 AL, DELTA! , SN,NR S K I P , L AP , NC ROS S , I CROSS ( 2 , 20 ), NCHP , YL A BE L ( 21 , I wlNOOw 
l,Fl,E2,ITYP£SP, RCGN,NESKIP , lEF, LAGi, LAG2 
C0MM0N/8LK2/irH( 14 I , CHSUM ( 14 » , N )FF( 14 ) ,CHSUMSU( 14) , S IGNA ( 14 ) ,RMS ( 1 
14), MEAN! 14) ,SCALFAC( 14) , C HSUM 1 ( 1 4 ) , TR AC K, ! 14) , ICHANlW) 

COMMON/ BLK3/NPT ,TMAX , NPT02 , NS PCT , UElF ,N o4,NPT012B 
1, INZERO,NREAD 

CJMM0N/SLK5/PCTC,NBINS, OMAXI 1 4 ) , UMIN 1 1 4 ) , UB 1 N ( 14 ) , B IN S I 1 00 , 14 ) 
COMMON/BLK6/IS AVE64, IR I , I X PLOT , I U ATA , 1 Z , 1 SPEC! 

C0MM0N/BLK8/I AGTGSP! 14) , 1 AUTOCO! 1 4) , 1 OR S P !20 ) ,I CRCOR !20),ITRA!20), 
1ICOH!20) 

CCMMON/BLKS/NE I L TP , F RE OF ! GO ),WGt)TF!50) 

C 

C COMPUTE AND PRINT SPECTRAL FILTER 

IFINFILTP.LE.O) GO TO 1 

CALL SPL INE !FREQF,WGHTF, NF ILT P, NS PCT , CM A 1 N ! I R I ) , UE L F ) 

NREC = NBLK')‘NChP«- 1 

CALL WRITES !8 , CMAIN! IR I ) ,NSPCT,NkEC ) 

DO 2 I=1,MSPCT 

2 CMAIN! IXPLOT41-1 |=!1-1)*UELF 
Il=IXPLOT-l $ I2=1RI-1 

PRINT 900, ! I , CMAIN! II P I ) , CMA IN! 12 + 1 ) , 1= 1,NSPCT ) 

900 FORMAT! //«! SPECTRAL FILTLR WEIGHTING F J NC TI ON*/ / 5X* 1* 3 A=!‘FhEOUENC Y* 
16X*WE 1GHT*3 !7X* 1*3X*FREQUENCY*6X* WEIGHT * ) / ! I 6 , 2 E i 3. 8 , 1 6 , 2E 1 3. !> , 1 6 , 
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22E13.5,I6,2E13.5)) 

1 CONTINUE 

CALL AUT0(CM.AIN ( Ul ,CM.AIN{ ISPECT ) ,CMAIN ( I «I) ,CMAIN( IXPLOT I I 

CALL AUTO EUNCTION ROUTINE WITH COMPUTED BLOCK ADDRESSES 

IF( NB INS .GT .0 > CALL NOTHAL 
RETURN 
END 

SUBROUTINE AUT 0 ( Z t SP E CT ,P. 1 ,XPL0T > 

DIMENSION SPECT ( II ,XPLQT( 1 » ,RI( I J ,Z( 1 ) 

COMPLEX Z, SPECT 

COMMON/BLK1/STARTT,ITFMT,NBLK, I PO W2 , NCR , NPR I N T , I PL OT A , IPLJTC , OFF SC 
iALf DELTA T, SNiNRSKl P, LAP, NCROSS, ICROSS (2 ,20) , NCHP , YLABE L ( 2 I , I WINDOW 
l,Fl,F2, ITYPESP,WCUN,NFSKIP,IFF,LAG1,LAG2 
C0MMQN/BLK2/ICri( 14) ,CHSUM( 14) ,NUFF{ 14 ) , CHSUMSQ) 14 ) , S 1 GO A ( 1 4 ) , KM S ( 1 
14) , MEAN! 14 ) ,SC ALFACI 14 ) ,C HSUH 1 ( 14 ) , TH AC K ( 14) , I CHAW I 14), I FILTER 114) 
C0MMGN/8LK3/NPT ,TMAX ,NPT02 ,NSPCT , OE LF , Nt 4 , NPT 0 1 28 
1 , INZERO, NR FAD 

C0MM0N/BLK8/I AUTOSPl 14 ) , I AUTOCOl 14),ICRSP(20),I CP.COR (2u) , ITRA120 ) , 
1IC0H120) 

CGMM0N/BLK9/NE ILTP 
COMMON/ BL-K 10/ NR CRD Z 

DIMENSION PLABELl6),PPLOTI27),BANDl2n , IDENIS) 

REAL MEAN 
COMPLEX ZZ 

PRINT 912, (TR ACKl J) , NOFE 1 J ), J = 1 ,NCH) 

912 FORMAT!//* NO. CF OFFSCALE VALUES FUR EACH CH ANNEL */ 1 1 X , A 1 0 , r 10 ) ) 

COMPUTE MEAN AND VARIANCE OF EACH CHANNEL 

00 36 IKK=1,NCHP 
IK= ICHAN ( IKK) 

MEAN! IK)=CHSUM1(Tk)/(NOLK*NREAO) 

SECMOM = CHSUMSC( IK )/ ( NBLK*NREAD) 

SIGMAS0=ISECMGM-MEAN( IK)**2 ) 

IF(SIGMASQ) 38,37,37 
38 PRINT 913, IK 

913 FORMAT!* NEGATIVE SIGMASQ FOR CH.,*I5) 

S IGMA! IK ) = 0. 

GO TO 36 

37 SIGMA! IK )=SQRT( SIGMASQ) 

36 CONTINUE 

START OF LCOP FCR COMPUTING AUTO SPECTRA AND CURRELATIUNS 

DO 50 JJ=1,NCHP 
J=ICHAN!JJ) 

DO 53 I = 1,NSPCT 
53 SPECTII )=0. 

C AVERAGE NBLK SPECTRA FOR ONE C'-'ANNEL 
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00 5^ IBLK=1,N|BLK 
IJ'=( TBLK-1 )*NCHP+ JJ 
CALL REAUMS (8, Z ,NPT, IJ» 

GO TO (80,80, 81 ),ITYPESP 
80 00 91 I=1,NSPCT 

91 SPECTd ) =SPFCT ( 1) t(P.EAL(Z ( I n'i='*-ZtAIMAG( Z( 1 ) )«*2» 

GO TO 84 

61 on 82 T=l,NSPrT 
ZZ=Z ( ! I 

AMP=f.A8S 1 ZZ ) i ARG=0. t IE ( A,v,P .NE . 0. J ARG = AT AN 2 ( A I MAG ( Z Z J , RE AL ( Z Z » ) 

82 SPECT( I ) =SPEOr ( I ) tCMRLX ( ARP, ARGI 
54 CONTINUE 

COMPUTE Cr-KRECTEG AUTO SPECTRUM ANO CALCULATEO VARIANCE 
SUM=Oo 

GO TO ( 8.1, 83,84) ,TTYP£SP 

84 C0N = 2o / ( NBLK^^NR E AO) 
on 85 I=1,MSPCT 
AMP = P,EAL (SPEC! ( I ) )<'Cn.N 
APG= AI MAG( SPECT ( I ) ) « 5 7 . 29 5 7 79 ^/NRLK 
SPECT( I ) =CMPLX ( AMP,ARG) 

85 CONTINUE 
RMS ( J ) =0. 

GO TO 86 

83 CfiN = DELT AT*0ELTAT/(6.2a3185308<n-.C0;'i’»fi8LR) 

00 59 I=1,NSPCT 
SPECT{ I ) =SPECT (I )*CON 

59 SUM=SUM*-SPECT ( I ) 

RMS(J) = SqKT(SUM*t:ELE *12.56637062) 

86 CONTINUE 

C APPLY spectral HLTER 

IE(NEILTP.LE.O.CRoIFILTEK( J),EO.L) GO TO 58 
CALL REA(JMS(8,RI ,N SP C T , NB L K*NCHP » 1) 

00 112 I=1,NSPCT 
112 SPECT( I ) = SPEC T ( I )*RI ( 1 ) 

C 

c PRINT AUTO SPECTRUM ANO SET Ul^ PLOTTING ARRAYS 

C 

58 PRINT 1001, J , TRACK! J) , ME AN( J ) ,S I GMA( J) ,RMS ( J) 
iOCl FORMAT ( ///*1CHANNEL* ! 3 , 3XA10, 3X*M E AN=* E 1 2 . 5 , 3X*SIGMA=* 

1E12.5,3X*PMS=*E12.5) 

1002 Ff)PMAT(///^ AVERAGE UF*I3,* T RAN SFORM 5* / / 5X * I *3 X*FR EUU ENC Y*oX *POWE 
1R*3 ( 8X*I *3X*FR ECUENC Y*6X*P0WER*) ) 

IF( ITYPESP.EO. 3 ) GO TO B7 
OU 60 1= 1,NSPCT 
XPLOTI I ) =( 1-1 ) *CELF 
AMP=SPECT (I )*2. 

SPECTI I ) =R I ( I ) =AMP 

60 CONTINUE 

911 FORMAT! 4( I 6,2t 13.5) ) 
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CALL v.'RITMS(9,RI,NSPCT, JJ ) 

IFC IAUTDSP( J).LE.O) GO TO 100 
GO TO (66,67,67 ), ITYPESP 

66 ENCOUL ( 60, 901 , PLA8EL ) L)E LF , TP ACK ( J ) 

901 F0RMAT(«AUT0 POkER SPECTRUM ( SANDw lUT H= * E9. * I *9X , A 10 » 

NPLA8EL=A1 

DO 69 I=1,NSPCT 
69 RI( I ) = RI ( T IX'DELF 

PRINT lOOS, (PLABEL( I » ,I=1,5» 

GO TO 68 

67 ENCOOE(bO,902,PLARELI TRACKUI 

902 FORMAT(«AMTn POWER SPECTRAL DENS IT¥*3 a, A i0,20X) 

NPLABEL=27 

PRINT 1006, (PLABEL( I) ,1=1,31 
1005 FORMAT! //20X7H * « ^ ,5A10,6H * * 
lOOo FORMAT ( //20X7H * >ir ♦ ,3A10,6H ♦ « 

GO TO 68 

87 00 88 I=1,NSPCT 
XPLOT! ri=(I-ll*OELF 

88 RM n = SP ECT ( I I 
ENCODE(60,903,PLABEL) TRACK(J» 

903 FORMAT!^ AMPL IT LDE SPE CTRUM AND PHASE * 5 X , A 1 0 , 16X ) 

PRINT 1009, PLASELd ) ,PLAeEL(2» 

1009 FQRKAT(//20X7H «■!■■(= ,2A10,6H « * 

NPLABEL= 18 

PRINT lOOd NBLK. 

1CC7 FURMAK///* AVERAGE 0F*1^,* TRANSF0RaS*//5X<=I *JX«FRE0UENCY’‘'4Xt‘AMPL 
1 ITUDE’t'4X*PHASe*2 ( 9X * I ’('SX f-FREuUENCY^^X* AMPL 1 TUUE ♦4X«RH A S E* J ) 

PRINT 1008, ( 1 ,XPLOT( n ,SPECT(I» ,I = 1,NPRINT) 

10C8 FORMAT ( 1 6 , 2 El 3 . 5 , F8, c , 1 1C , 2E I 3. 5 , F 8, 2 , 1 10 , 2 1 1 3. 5 , F b, 2 ) 

GO TO 89 
6C CONTINUE 

PRINT 1002, NBLK 

PRINT 911, ( I ,XPLCT( I ) , RI ( I ), 1 = 1 , NPRINT ) 

89 CONTINUE 
NRCRD7=NRCRD7f 1 

WKITE(7) •PLABEL,NSPCT,(XPLOT(I),Rl(H,l = i,NSPCr) 

PRINT 1010, NRCR07,PLAbEL 

1010 F0RMAT(//21H ♦ » * * « RECORD NO. ,15, ON TAPE7 CONTAINS *6A10,1G 
IH 4 « * * *1 

C PLOT FANFULO PLOTS AND COMPUTE 1/3 C.8. SPECTRUM 

GO TO ( 7C, 70,62 ,63,62,63 », I PLUTA 

62 ILOG=0 $ GO TO El 

63 IL0G=2 

6 1 CALL PLOTNBIYL ABEL ,TRACK( J I , 10,XPL3T,RI , NSPv-T , I LOG , F 1 , E 2 , P L AB EL , 
1NPLA8LL, IFE,1 ) 

IF( ITYPES°.LT.3) CO TO 70 
ILOG=C 

DO 90 I=1,NSPCT 

90 RI ( I l=AIMAG(SPECT(I ) ) 

ENCODE (50,904, PLAliEL I 
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9C4 FORMAT( <=PH4Se ANGLE, D£ GR E E S ’i'iOX ) 

NPLABEL=20 

CALL PLaTNi6(YLABEL,TPACK( J) , 1 0, X PlQT , R I ,\SPCT , I LOG , F 1 , F ^ , PL Abt L , 
INPLABEL , 0,0) 

GO TO 100 
TO CONTINUE 

IF( ITYPESP.E0.3 ) GO TO ICC 

COMPUTE 1/3 OCTAVE BAND SPECTRA 

GO TO (100, 6S>, ICC, 100, 65,65) .fPLOTA 
65 CALL HANDS ( OE L F ,NSPC T , S P bCT , P PLOT , BAND , i\j pp, 1 fc RR , I T V P t S P ) 
IF(NPP.LE.O) go to 100 

CALL PLOTBI YL ABEL, TRACK! J ) , 10 , 8 A NO, PP LO T , NPP , PL AbE L , NP L AcEL , IFF ) 
64 CONTINUE 

GO TO (72,71) , ITYPESP 

71 ENCODE (60, 73, PL ABEL) TRACK(J) 

73 FORMAT(*l/3 OB ALTO POWER SPECTRAL DE NS I T Y«6 X , A 1 , 1 OX ) 

GO TO 74 

72 ENCUDE(60,75,PLABEL) TRACK(J) 

75 FQRMAT(=S'l/3 OB AUTO POWER SP EC T RUMVAX , A 1 0 , 2C X ) 

74 CONTINUE 
NRrR07 = NRCRU74^ I 

WRITE(7) PLABEL ,NPP, ( BAND ( I ) , PPLOT( I ) ,I = i ,NPP ) 

PRINT 1010, NRCRDZ.PLABEL 

100 IF( I AUTOCU( J) .LE.G) GO TO 50 

compute AUTOCORRELATION 

00 101 I=2,NSPCT 
Zn)=SPECT(l) 

101 Z(H-NSPCT) = SPFCT(IjSPCT-I»2) 

Z(1 ) =Z(NSPCT«-1 ) =SPECT ( 1 ) 

IF( ITYPESP.GT. 1 ) GO TO 113 
DO 114 I=1,NPT 

114 Z(I )=Z( I )/DELF 
113 CONTINUE 

CALL FOURT( Z, MPT, 1,-1, 0,S°ECT) 

C0N=3. 14159265 3 58SB/TMAX 
IFIINZERO) 1014,1014,1011 
1014 DO 110 I=1,NSPCT 
110 SPECUl )=Zi D^-CCN 
GO TO 10 12 

1011 DO 1013 1=1,NSPCT 

1013 SPECK I ) =Z { I )*CCN*NPT/ ( NPT-I-I +2. ) 

1012 CONTINUE 

PRINT AUTOCORRELATION 

D(3 102 1 = 1, NPT 
IM1=I-1 


101 



o o o 


APPENDIX H 


RI( I » = SPECT(1 ) 

102 XPL0T( 1 )=IM1 

PRINT 1003, TPACK( J) , (XPLOK I I ,K I (I ) , 1 = 1 ,NSPCT I 

1003 FORMATC //^lALlTO CORRELATION , *A 1 J// 8 ( 7X I >:‘oX>»R ) / ( 2X , F ‘j .0 , E 1 1 . 3 , 

1F5.0,E11.3, F3.0 ,Elio3,F'j.C, Ello3,F5.C ,E 11.3,E3. C,Elli,3,Eb,0,Eli.3, 
2F5.0,E11.3n 

ENCODE (60, 10C4, FLABELl TRAf.K(J) 

1004 FORMAT(«AUTO CORRELATION *,A10, 32X) 

NRCR07 = NRCRD7<-1 

WRITE(7> PLABEL ,NSPCT, ( XPLOK n,RI( I ) ,I = 1,NSPCT I 
PRINT 1010, NPCRL7,PLABEL 
IF(LAG1.E0.0.ANC.LAG2.E0.0) go to 50 

SET UP PLOT ARRAYS AND PLOT AUTGCURRE LAT 1 ON 

L1 = LAG1 $ IF(Ll.LT.OI L 1 = 0 $ L2=LAG2 
IF(LAG2.GT.0) GO TO 109 
L1=-LAG2 i L2=-LAG1 
109 I1 = LK1 i I2 = L2+1 

106 IF( I1.GE.I2I GO TO 50 
NPLT = I2-I1+^1 

OQ 107 I=1,NPLT 
11=1+11-1 
XPLOT( I »=I r-1 

107 RI( I l=SPECT(II I 

CALL ASCAL E(R 1 , LO.,NPLT , I, 10. » 

IF( IFF.EO.O ) GO TO 50 
IF(NPLT-25j» 201,201,202 

201 K=1 $ GO TO 103 

202 K=NPLT/256 

103 ENC0OE(50, 920, IDEM TKACK( J) 

920 FORMAT(*AUTOCCFRELAT I0N*3X,A10,22X» 

PN=RI ( NPLT+ 1» 

PX=RI (NPLT+2) «10.+PN 

CALL FANFOLD (P I ,NPLT ,K, 1 , NPLT , I OE N , I H* , 1 . , PX , PN , Y LAB EL , 2 , UO , 0 , o. , 
IIH I 

50 CONTINUE 
RETURN 
END 

SUBROUTINE NORMAL 

C0MM0N/BLK1/STARTT,ITFMT,N8LK, I F Ow2 , NC H , N PR I NT , I PLOT A, IPLOTC , OF ESC 
1 AL, DELTA T, SN, NR SKIP, LAP, NCROSS, I CROSS 12, 20 » , N CHP , YL ABE L ( 2 ) , IWiNOOw 
1,F1,F2, ITYPESP 

C0MM0N/BLK2/ICH( 14) , CHSUM( 14) ,NOFF( 14 ) , CHSUMSQ( 1 4 ) , S I GM A ( 1 4 ) , RM S (1 
14),MEAN( 14) ,SCALFAC( 14 ) ,C HSUM 1 ( 14 ) , TRAC K ( 14) , ICHAN( l4) 
.C0MMCN/BLK3/NPT ,T MAX , NP TU2 , NSPC T , UELF , N64 , NPTO 1 26 
1 , INZERO, NRE AO, NFTOT 

COMMOM/8LK5/PC TC,NBINS,CMAX( 14) , DMI N ( i 4 ) , UB I N ( 14 ) , B IiM S ( 1 00 , i 4 ) 
1,CHIS0C 
REAL MU, MEAN 
DIMENSION YL(3) ,ICEN(5) 
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nATA IDEN/ICHCQUNTS ,4*1.0H / 

NbINS>i3 = N8 INS-3 
00 1 T=1,MCHP 
J=ICHAN(1I 
MU=MfcAN( J) 

ENCODE(30,9Ci;, YL) TRACMJ) 

902 EORMATI 3X*HIST0GRAM EOR *A10,3X) 

PN=PX=d, 

CALL FANEn|_o( HI\S( 1, I I , MB INS , If 1 ,lCOi I JEN, Ih*, 1. ,PN ,PX , YL, 3, 
llOOfCfXAPPAYfXLABELI 
SIG=STGMAU) 

EACT08=1,/ (2o80662B2 
OFL = DBIN { I) 

0MAX=0. 

SUH=C. 

lElOMINl J>,LT,-20. ) TO 9 

P,N= PFUN ( -20. , OM { ,J » , MU , S 1 G ) «EAC TOR 

CO TO 7 

6 PN=0. 

7 CONTINUE 

PRINT 905, (B!NS(K,I l,K=l,NBINSI 
'9C5 FORMAT! /« B I N / ( iOF 13 . 0 ) ) 

CHI S0=0. 
on 2 K=1,NBINS 
A = Kt<DEL«-UMIN( J ) i 

SUM= B!MS(K,n . 1 

P=SUM/NPT0T j 

PN= PFUN(A-U£L,A,mk,S IG)*FACTOR 
IF(Pi:.EO.d) GO TO 2 
CHI SO=NPTnr*( P-FNt RNfCHISO 
2 CONTINUE 

ALPHA=l. -RCTC/ICC, 

PRINT 900, N8INSM3,CHIS0, ALPHA, CHISQC 
90C FORMAT (/// 12X 13F « ■c< GOrONESS OF FIT TEST * * «'//<' UEGKElS UF FR 

lEEOQM =i«I5,1CX’!'CHI-S0UARL =#F1C,3//^' AT THE SIGNIFICANCE LEVEL OF* 
2F10.3, ICX^THE CRITICAL VALUE OF CHI-SOUARE IS=>F10.3) 

1 CONTINUE 
RETURN 
ET'O 

FUNCTION PFUNI A,P,MU,S!G) 

kfal mu 
PFUN=0. 

0X= ( fl-A ) /50. 

X = A-[)X/2. 

00 1 1=1,50 
X=XfOX 

PFUN=PFUNfI)X«lXP(-(X-MU l*>!=2/I 2.*S IG^'SIG I ) 

1 CONTINUE 
RETURN 
END 

SUBROUT I NE PLUTP (VLABEL , F R AMEL , NF , BAND , P FL OT , NPP , P L ABE L , NP , I FF ) 
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OlMENSl uN 3AND( 1 ) , Pt>LOT ( 1 1 , YLAHtL <2 » , FKAMEK 1) ,PLABELl5),FFlU(5) 
1,I0EN(6) ,BCF(24) 

OATA HuF/50*fO3*i tiO* vIOO* fl23»tl30« f20D« t230*t313*f^00*y 300 • y 03 w • t 
IdOC* f 1000* f L2 30 at lOOOa t 2000* f 2300«t 313Co t^GOO© f dOOO© yBOuG© y lOOGC© y 
220000./ 

plot POWEk fop EANU CEMEK FREUUENCIES 30-2CK HZ 

ACOl DO 2001 1 = 1 yN'PP 

2001 BAND(I)=I+A 
PMAX = PPLOT( 1 I 
00 2002 I=2yNPP 

2002 PMAX = AMAXl(PPLQTn lyPMAX) 

NMAX=IFIX(PMAX) 

IF1NMAX.lt. PMAX) NMAX=riNlAX«-l 
PMAX=NMAX 

PMIN=PMAX-3, 

I2=NPP 

on 2004 I=lyI2 

IF( PPLOT ( I) .LT. FhINI P P LOT ( I » = PM I N 
2004 CONTINUE 

PPLQT (NPP«-n = PM IN $ PPLOT (NPPV2) =.3 
BAND(NPP*1 ) =1. t HANOI NPP ♦■2 ) = 3. 

IF( IFF.FQ.O) RETUPN 

F.FIOm = YLABEL ( 1 ) i FF I 0 ( 2 > =YLAb EL ( 2 ) 

NWOPOS = ( NF >9) / 10 
DO 1' 1 = 1 yNwOPOS 
1 FrID(l<-2 ) = FRAMEL ( n 

ENC00E(34,9C3, lOENJ PLABEL 

905 FORMATIOLOG *5A10I ... 

CALL FANFOLOIPPLCTyNPPy ly lyNPPylOENy lH»y l.yPMAXyPMlNyFFlUyNWORUSft 
1 y 120y lyBCF , ICHFREwUENCY I 
RETURN 
ENO 

SUBROUTI NF BANOSIOELFy NSPCTy SPEC T, PPLOT y BAND ,NPP y I E KK , I PGWPL T ) 
COMPLEX SPECTIll 

DIME NSr ON' FPLOT (1 ) , IND(46,2»yFNC( 4fa» ,FNL (4 7) 

OATA FNL / I. 1220 , l.4i2 5 , 1. 7 783 y 2. 2 387 , 2 . 8 1 84 y 3 . 5 48 1 y 4. <y8o 3 y 5 .o23 4 , 

1 7,C793y 8.9 123y 1 1.220 , 14. 123y 17. 783y22© 3 8 7,28. 184 , 33. 48 1 , '♦4 . 8t>8 , 

2 50.234, 7 0. 7 S3, 89. 123, 11 2. 20, 141.^5, 17 7. 3 223 . 8 7 , 28 1 . 8 , 334 . 81 , . 

3446. 6 8,5 62© 34, 7 C 7© 95 , 89 1 . 2 3 , 1 12 2. 0, 1 4 1 2 . 5 , 1 7 7 3. 3 , 22 38. 7 , 28 18 . 4, 

43 348. 1, 4 466. 8, '3623. 4 , 7 C 79 . 3 , 8;1 2. 5, 1 i 22 0 . , 14 1 23 . , 1 7 7 83 . , 223 a 7 . , 

5 28184. , 3 5481. ,44668. / 

DATA FNC/1 .23 89 ,1. 3849, 1 . 9953 , 2. 5 1 1 9 , 3. 1623, 3. 98 1 1 , 5 . 0 1 1 9 , 6. 3096 , 

1 7.9466, IG. , U. 3 69, 13 . 849 , 1 9. 9 53 , 2 3. 1 19 , 3 1 . 02 3 , 39 . 3 1 1 , 3 C . 1 1 9 , 63. 096 
2, 79.433, 100. , 125.09, 156.49, 199.53,251. 19 , 316. 23, 398. 1 1 , 5G1 .19,, 

3 630. 96, 7 94. 33 , ICCC . , i 2 3 8. 9 , 13 84. 9 , 1 99 5. 3 , 23 1 1 . 9 , 3 16 2. 3 , 3,98 1 . 1 , 

4 5011.9, 6 309.6, 7 943.3, 10000. ,12589., 13 84 9 ., 1995 3. , 2 3 1 19 ., 31 62 3. , 
339811. / 

DIMENSION PSr)(46),PSDPHZ(46),8AN0(27) 

Dli 3 1 = 1,? 7 
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3 B;iNl) ( I )= I 
IER« = 0 
DO I 1=1, ^6 
BMD = FNL( H-1 )-FNL( 1 ) 

I 1= I 

IF( BNU.GE.OELF ) GC Kl ?. 

1 CQNTTJUE 
IERP=1 
RETURN 

I DO 62 1=1,40 

PSDPHZi I l = PSO( I )=G, . 

62 INO( 1, 1 ) =IND( I , 2 (=0 
I 1= t 1 

l)0 63 ! = 2,N5PCT 
F=J I-D’i'DELF 
68 IFIII— 46) 64,64,65 

65 I I = 46 

GO TO 66 

o4 IF(F-FNL ( I I +1 ) ) 66,67,6? 

66 IF ( IMO( I I , 2 )o EQ.G) lr;0(ll,l) = l 
IND( II,2)=IND( I I,2M-l 

GO TO 63 

67 ii=in-i 
GU TO 68 

63 CONTINUE 
NPP = 0 

■ iJO 76 1=11,46 

IF( INCH I ,2 I » 76,76,7 J 
70 ISTRT=INO( I ,1 ) 

NPP=NPP+1 

- CALL RNOSlIM (SPECT ( ISTRT I , IND( 1,2 I ,PSUl I n 
PSOn )=PSD( I) ^DELF, 

PSOPUZ ( I ) = PSG ( I )/ IFNL ( I n )-FNL( I ) ) 

76 CONTINUE 

NPP=NPP-( 17-1 II 
IF(NPP.LT.O) NPF=0 
IF( NPP.GT.27) NPP=27 

PRINT 80, (FNCI II ,PSl)( I l,PSUPHZ( I I, I=Il,4o) 

80 FORMAT!//* 1/3 OCTAVE B 7 NO* 5X* P O'wE F *6 X* POUF K S P ECT,'. A L * / * CENTER 
1 fRE0UENCY*4X*SFECTPUM*7x*0EN5lTY*/(Fl2.C,E19,4,E15o4) I 
IFINPP.EO.O) GO TO 3004 
GO TO (3001, 30021, IPOWPLT 
3001 00 30C3 1=1, NPP 

PPLOTCI I =-100. t IF! PSD ( 1 ♦•16) .GT .0. » PF LOT! I I = ALOGiC ( P S D ! I f 16 I I 

3003 CONTINUE 
GO TO 3004 

3C02 00 3005 1 = 1, NPP 
PPLOTd I =-100. 

IF! PSUPHZ! I *16 » .GT.O. I PPLOT ! I I = A LOG 10 ! PSUPHZ ! I *16 I I 
3005 CONTINUE 

3004 RETURN 
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END 

SUBROUTINE BNU SL'I^I S , Njyi , P SO » 

COMPLEX S(l) 

PSD=0. 

IFINUMoGTe 1 ) GO Tt! 1 
PSD = PEAL (S( IJ ) 

RETURN 

1 PSD=REAL (SI 1) »<-REAL(S(NUM) ) 

IF(NUM.EQ.2I RETURN 
NUMM1=NUM-1 

00 2 I=2,NUMM1 

2 PSD=PSDvREALtS( 1) ) 

RETURN 

END 

SUBROUTINE SPLINE ( X , Y , KNT , KNTOUT ,Y0UT , OX ) 
DIMENSION C (^,4 > , A(^ I , IP IV( 4) ,X( I ) , Y( U , YOUTI l» 
XM = 2,fX( I ) - X(2I 
YM = 2.=fY( 1) - Y(2I 
XN = XII) 

YN = YU 1 
XU = XI 2) 

YO = Y(2 » 

XP = X(-3) 

YP - Y(3 ) 

SLN = ( YN-YM) / ( XN-XM ) 

SLO = (YO-YNI/ (XO-XN) 

SLP = ( YP-YO W ( XP-XC ) 

RKQUT = KNTOUT 
IYNX=X(1I/0X 

IF( IYNX*0X.LT.X ( 1) ) I YNX=IYNXU 
VAR= I YNX*OX 
DO 1 I=l,IYNX 
1 YOUTH 1 = Y( I ) 

I YNX= IYNX+ 1 

LIM = KNT V I 

DO 7000 N=3,LIM 

IF (SLN ,NE, SLC .OR. SLO .NE. SLP) GO TO 30u0 
C LINEAR 

A(AI = 0. 

A(3) = 0. 

A(2) = SLO 
All) = YN - SLO«XN 
GO TO 60GC 
3000 CONTINUE 

Ad I = (SLO - SLN)/(X0 - XM) 

A(2) = (SLP - SLO)/(XP - XN) 



IF (A(l) 

,NE. A (2 ) 1 GO TO 5000 

c 

PARABOLIC 



C( 1,1) = 

1. 


C(2,l) = 

1. 


C(3,l) = 

1. 
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CU,2) = XN 
C(2i?l = XO 
C(3,2) = XP 
C(l,3» = XN*Xi\ 
c(2,3) = xn=!=xn 
G(3,3» = XP*XP 
AU ) = YN 
A(2I = YO 
A(3) = YP 

call SIMLQ (C, 3t A, 1 >, DLT, IPIV, I SC » 
,AK) = 0. 
cn TO, 60 GO 
5 000 CfJNTU.'UE 
C CUb I C 


C(l,l) 

- 

1 . 

C(2,l) 

= 

1. 

C(3, 1) 

= 

0 • 

C(4,l) 

= 

n. 

C( 1 ,21 

= 

XN 

C( 2 ,2) 

= 

xn 

C(3,2) 

= 

lo 

C(4,2) 

- 

I . 

C(l,3) 


XN*XN 

C ( 2 , 3 ) 

= 

xn^xu 

C(3,3) 

= 

P.^XN 

C(4,3) 

= 

2. «XO 

C(l,4) 

= 

XN«C (1,31 

C(2,4) 

= 

xn*C( 2,3) 

C(3,4) 

= 

3«*C( 1,31 

C(4,4) 

= 

3.*C(2,3) 

A( 1 ) = 

YM 

A(2) = 

YO 


A(3I = TAM{,5-!‘{ATAM(SLN) + ATAN(SLU)n 
A(A) = TAN(.5''( ATA>USi-OI <■ ATAN(SLP)J) 

CALL SIMEO (C, 6f A, I, iU T , TPIV, 4, ISC) 
6000 CONTINUE 

SUM = A ( 1) 

VARP = 1. 

00 6100 K-^2,4 

VARP = VARP*VAR 

SUM = SUM A(K)*VAEP 
6100 CONTINUE 

YOUTIIYNX) = SUM 
VAR = ELOAT( I YN X ) *UX 

1 YN X = I YN X «• 1 

IF (VAP .LE. X(N-in GO TO 6000 

XM = XN 

y.M = YN 

XN = XO 

YN = YO 

XO = XP 
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YU = YP 

IF (N .LT. KNT) GO tq <,500 
IF (N .EQ. KNT«-1) GU TO 7000 
XP = 2.«X(KNTI - X(KGT-l) 

YP = 2.*Y(KNT> - Y(KMT-l) 

GO TO 66 QG 
65C0 CGNTIMUE 

XP = X(N+1 » 

Y P = Y ( N 1 1 ) 

6600 CONTINUE 
SLN = SLO 
SLC = SLP 

SLP ,= (YP - YCl ) /(XP - XO) 

7000 continue 

IF( I YNX.GT.KNTQLT) RETURN 
DO 2 !=1 YNX.KNTCUT 
2 YOUTH ) = Y(KNT) 

RETURN 

END 

OVFRLAY( PATS,4,GI 
PROGRAM CPOSSRP 
COMMON CMAIN( H 

COMMON/HLK 1 /STAR TT, I TFM T , N8 LK, I P0W2 , NC-t , NPRINT , I PLOT A, I PLUTC.OFFSC 
lAL.UELTAT, SN,NRSK:P,L AP.NCRUSS, ICR0SS(2 , 20 » , NCHP , YL A 8E L ( 2 1 t IwINUO;/ 
lfFliF2, (TYPES P, RCuN, NFSKIF , IFF , LAGl , LAG2 
C0MM0N/BLK2/ICH(14l,CHSUM( 14 I , N OFF (14 ) , CM SUM Sw ( !*♦) ,SIGMA(14) ,RMS(1 
14),MEAN( 14) ,SCALFAC( 14> ,CHSUM1( 14) .TRACK ( 14) ,ICHAN(14) 

C OMMON/R LK 3 /N FT ,TMAX ,NPT02, NS PCT ,OELF ,N 64, NPTOl 2b 
l.INZERG.NREAD 

CGMM0N/BLK6/I SAVE64, (RI , I XPLOT , lUATA, U, ISPECT 

CQMM0N/6LR8/1 AUTOSP( 14) , 1 AUTOCO( 14) , IC R S P (20 ) , I C RCOR ( 20 ) , I T RA ( 20 ) , 
1ICOH(20) 

CQMM0N/BLK9/NF ILTP 
ccmmon/rlkio/nrcro? 

CALL CROSS ( CMAINH/I ,CMA IN( ISPECT ),CMAIN (IRI ) ,CMAIN( IXPLOT ) ) 

RETURN 

END 

SUBROUT! NE CROSS H.SPECT ,R 1 .XFLOT ) 

COMPLEX SPE CT ( 1 I ,Z( 1 I 

DIMENSION P 1 U ) iXPLOM 1 ) ,ERAMEL( 3 ) , PLABEL (6 ) 

C0MM0N/BLK1/STAPTT,ITFMT,NBLK, I PO W2 , NC H , WPR 1 N T , 1 PL OTA, IPLOTC ,OFESC 
1AL,0ELTAT,SN,NR SK I P , L AP , NCROS S , I CROS S { 2 , 20 ) , NC HP , Y LABE L ( 2 ) , I^ INOOrt 
1 ,F1 ,F2, I TYPESP,>CON,NFSk IP , IFF,LA(,1,LAG2 
C0MM0N/B.LK2/ICH ( 141 , CHSUM( 14) ,NOFF( 14 ) , CHSUMiO( 14) , S I GMA ( 14 ) , RMS ( 1 
14),MEAN(14),SCALFAC(14) , C HSUM 1 ( 1 4 ) , TR AC K ( 14) , ICHAN( 1 4 ) , I F I L TE R ( 14 ) 
C0MMGN/BLK3/NPT , TMAX , MP T 02 , NSPCT , UELf ,No4,NPT012B 
1 , INZEPO 

COMMON/BLKb/I AUT0SP( 14) , I AUTOCOI i 4) , 1 C KS P 1 2 *: ) , I CKCCiR(2L ) ,-ITRA (20 ) , 
1ICOH(20) 

C0MM0N/6LK9/NF 1 LTP 
COMMON/BIK 10/NR CRD7 
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DIMENSION I DEN (5), BAND (27) tPPL0T(27) 

DATA RAD/57.2957795/ 

START OF LfOP FCR COMPUTING CROSS FUNCTIONS 
DO 70 !CR=1,NCRCSS 

M= ICROSS ( 1 f I CR ) $ K2= ICRCSS(2,ICR ) 

IF( lAUTOS” ( K1 ) . EC.O.OK. r AUTUSP(K2 ).eo.O ) GO TU 75 
00 77 I=1,NCHP 
IFIKl.EO.ICHARM I )) J 1= I 
IF(K2.F.Q,ICHAM I 1 1 J2=I 
77 CONTINUE 

DO 76 I = 1,NSPCT . 

76 SrECT(n=0. 

AVERAGE C^OSS SPECTRA FOP ONE PAIR 

DO 71 IBLK=1,N6LK 
IJ= Jl+( I BLK-1 ) *NCHP 
CALL KEADMS (8 , Z ,NPT, I J ) 

IJ= J2<-( I 3LK-1 ) *NCHP 

CALL RE ADMS (8 , Z (NSPC T«-1),NPT,IJ) 

DO 71 I=1,NSPCT 

71 SPECT( I )=SPECT( II«-Z( I l«CONJG(Z( I ♦■NSPC T) ) 

IF(NFILTP.LE.O) GO TO 111 

1F( IFILTER( Jl) .EO.O.ANO.IF ILTER( J2).ED.O) GO TJ 111 
CALL READMS(8,RI ,NSPCT,NBLK*NCHPfl) 

DO 112 I=1,NSPCT 
112 SPECK I ) =SPECT( I )<=RI ( I ) 

111 CONTINUE 

CALL WRITMS(9, SPECT,NPT ,NCHPfl ) 

COMPUTE CORRECTED CROSS SPECTRUM 

CON=DEl.TAT't<UEL TAT/ (6.2 83185308«WCGN«N3LK ) 

DO 115 I=1,MSPCT 
115 SPECTC I ) =SPECT( I )>:=CON 

ENCODE(23t900,FRAHEL ) TRACK (Kl) , TRACK (K2) 

900 F0RMA-^(A10,3H X ,A1C) 

IF(ICRSP(ICR).LE.C) GO TO 105 
GO TO (60,61) , ITYPESP 

60 CON=CON*DELF 
ENCODE(50,9Q1 .PLABEL I DELE 

901 FORMAr(*CRQSS PCWtR SPECTRUM (BANDWIDTH =*G9 . 2 , « )<■ , 6 X ) 
NPL ABEL=49 

PRINT 898, (PLABEK I I ,I=1,5),FRAMEL 
898 F0RMAT(>1‘1«5A1C,5X,3A1C) 

ENCODE (o 0,921, I DEN) TR ACK ( K 1 ) , TP ACK ( K2 ) 

921 FORMATIRCROSS POWER SPECTRUM *A 10, IX , A 10 , 7X ) 

GO TO 62 

61 ENCODE! 50,902 , PLABEL ) 
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902 FORMAK’i'C^OSS PCWER SPECTRAL JENS1TY*22X) 

NPLABEL=2a 

PRINT 899, (PLABELd ) ,1=1,3), FRAMEL 
899 FORMAT (*i=>3AlC, 5X,3A 10 ) 

ENCODE! 6 0,922, I DEN) TRACK(Kl) ,TPACMK2) 

922 F0RMAT(«CP0SS PCUER SPECTRAL DENSITY ^A10,1X,A10) 

62 IFINPRINT.GT.O) PRINT 897 

8 97 FORMAT ( / /4X * I A FRE JOE NC Y*9X*RE A L ♦ 1 C X* I M AO* 7 AMPL I TU DE <^7 X>^PHA S E=f 

1 ) 

SET UP PLOT ARRAYS AND PRINT CROSS SPECTRUM 

00 72 I=1,NSPCT 

F=( I-l )/TMAX 

XPLOT(I)=F 

AMP = CA8S (SPECT ( I ) ) 

IF(AMP) 79,80,79 
80 ARG=0. i GO TO 78 

79 ARG= AT AN 2 ( AIM AG ( SPECT ( I ) ) , REAL (SPECT ( I ) ) )*R4D 
78 CONTINUE 

IF( I-NPR.INT ) 73 , 73,200 
73 PRINT 915, I,F, SPECTd ) ,AMP,ARG 
915 FCRMATCI 5 ,5E14. 5 ) 

200 I F ( I PLOTC. EQ.O ) GO TO 201 

GO TO ( 201, 201, 202,202, 202, 202), IPLUTC 

201 RI ( I )=SPECT ( I ) 

RI( M-NSPCTf2)=AIMAG(SPECTd)) 

GO TO 72 

202 RI d ) = AMP 

RI( 1 +NSPCTf2)=ARG 
72 CONTINUE 

IF(IPL0TC-2) 301,301,302 

301 IDEN6-10H (REAL) 

I0EN7=10H ( IMAG) 

GOTO 303 

302 IDEN6=10H AMPLITUDE 
IDEN7=10H PHASE 

303 NRCR07 = NRCR07«- 1 
NSPCTP2=NSPCTf 2 

PRINT lOCC, NRCR07,IOEN, IUEN6 

1000 F0P,MAT(//21H ***** RECORD NO., 15, * UN TAPE? CONTAINS *oA10,10 
***** ) 

,(P,ITE( 7 ) IDEM, ICEN6, NSPCT , (XPLOT ( I ) , R I ( 1 ) , I = 1 , N SPC T ) 

NRCRD7=NRCRD7*-1 

PRINT 1000, NRCRC7, IDEM, IUEN7 

WRITE ( 7 ) I OEN, I 0EN7 , NS P CT , ( XPLOT ( I » , P 1 ( I *NSP CT P 2 ) , I = 1 , NS PC T ) 

PLOT FANFOLD PLOTS 

NN=(NPLA8EL-1 ) / 10*2 
NPLABEL=NW*10 
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T^^( I PLOTC. EQ.C ) r ,0 T ;1 ICb 
rj:1 TO ( 1 or , iC 7, ICS, lO'J, 1C V, 109) 1 PLOTC 
1C7 !I.OG = 0 i SO 10 1G6 
1C8 icnG=<^ 

106 ewcnuf.l lC,903 ,PlA8fL (NU) ) 

903 FORMftH * ( RE AU 

CALL PLQTOOl Yl. AEEL ,F R A M EL , 2 j , X PL OT , R 1 ,MSPCT, 1 L 0 G , E 1 , F 2 , P LA B t L , NPL A 
IBEL , IFF, 1 ) ■ 

F MC OOL ( 1 0 , 9 04 , P L A b EL' ( NO ) I 

904 FnP.MAT(^= ( INAG)-'^) 

CALL PLOTNB ( V L A EEL , F H AiB FL , 2 J, xPLuT ,P I ( \ S P CT f 3 ), NSPC T , I lOG , F 1 , F2 , 
1PLABE^,NPLA3EL , If-F,0 ) 

GO TO US 
1C9 ILOG=I PLQYC-3 
lu4 ENCOUEI 10,905, PL ABEi. ( Nw ) I 

905 FORMAT!* MAGNITULiE*! 

CALL FLuTNB (Yl APEL ,f ‘i AMEL , 2 J, XPl(IT,.-U ,.JSPCT , I LOG ,F 1 ,F2, PLABEL,NPLA 
13EL , ! FF , 1 ) 

ENCODE! 10,906,PLAf,El.!NW) ) 

906 FORMAT!* phaSL*) 

IF! IlOG.GT. 1) !LOG=ILOG-^ 

CALL PI OTNR !YL A EEL,F PA MEL, 23, XPLOT,R I !NSPCT*3 ) , NSPCT , I LOG, F 1 , F2 , 
1PLA3TL jNPLABEL , IFF,0 I 
105 CONTINUE 

IF I TCRCOR! ICR ) .LE.O) G(. TO 411 

COMPUTE CPOESCORRELAT ION 

DO 74 I=2,NSPCT 
Z!I)=SPECT!n 

74 ZINSPCTFI )=C0NJG(SPECT(NSPCT-I»-2) » 

Z( 1 ) = /. (NSPCTvl ) =;SPtCT I 1 ) 

IF! ITYPESP.GT. 1 ) GO TO il3 
on 114 1=1, NPT 
114 HI ) = 7I I )/DELF 
113 CONTINUE 

CALL F0UPT!7,NPT,1,-1,1,SPECT) 

PRINT 916, TRACK!K1) ,TRACK!K2) 

C0N=6.2 B31 H5308 /TMAa 
DO' 110 T = ), NPT 
no SPECTCI I =/. 1 I) *CCN 

SET UP PLOT ARRAYS, PRINT AND PLUT CURRELATICN 

IFllNZERO) 5CC,50C,5C1 
500 00 5C2 I=l,NSPC.T 
IM1=I-1 
IPN= 1 4-NSPCT 
RI ! I PN) = SPECT ! I ) 

XPLOT! I PN )= IMl 
IMN=I-NSPCT-1 
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R n n=SPECT ( I PINi ) 

502 XPLOTH ) = 

GO TO 303 

501 00 90 I=1,NSPCT 
I M I = T - 1: 

IPN= I *-NSPC.T 
I2M1 = 2*TM1 

RI { IPN'»=NPT*SPECTm / (NPT-I2MI) 

XPLOT( 1 PN) =IM1 
IMN= I-NSPCT-1 

I2MN = 2’«'IKM $ IFlI.tO.l) I2MN=I2MN*i 
R I ( I ) = IMPT=>SPEC T { IPM> / (NPTt-I2M.M» 

90 XPL0T( I I =rMM 

503 continue 

PRINT 917, (XPLCT(l) ,RI ( I ) ,1 = 1,NPT» 

917 F0RM6T(2X,e(F5.C,Ell.3) I 

916 FORKATI //*1CR0SSC0P.REL4TI ON, «A10,* X ♦ A li.// B I 7X * I 5X«RX 7* > ) 

ENCODE (60,923, PL ABEL ) TPACK(K1» , TRACk(K2) 

92-3 FORMAT! ^CKOSSCOPRELAT ION *A 1C , I X, A 10, 2CX I 
WRITE! 7 ) PLABEL ,NPT, (XPLDT ! II , R I ( 1 ) , 1 = 1 , NPT ) 

NR.CRn7 = NRCR 07 4-l 
PRINT 1000, NPCRrJ7,PLABEl 
IF! LACaoEO.Co ANC.LAG2.E0.0) GO TO 911 
I 1 = LAG1*-N9PCT + 1 t !2 = LAG2<-NSPCT+ 1 
206 IF! 1UGE.I2 ) GO TC 911 
NPLT=I2-I1*-1 

CALL ASC ALE!R 1( II ),10o, NPLT, I, 10. ) 

IF! ’FF.EO.G > GO TO 911 
IF(NPLT-256I 101,101,102 

101 K=1 $ GO TO 103 

102 h;=NPIT/256 

103 ENCODE! 50,920, IDEM) FRAMEL 

920 FORMAT!*CROSS CORREL AT I 0N*3X, 3A1 0 I 
PN=RI (NPLTf 1 1 ) 

PX=RT(NPL'"*-IIH )*10. fPN 

CALL FANFOLD! R I < I 1 1 ,NPLT ,K, I, NPLT ,lUhN,lH*,i.,PX,PN,YLABEL,2,12C,0 

1,0. ,1M ) 

GO TO 911 

75 PRINT 918, TR A C K ! K 1 1', T R AC K ( K2 I 

918 FORMAT!//* CROSS SPECTRA, EUR *A1C,* X *A10,* CANNOT BE COMPUTED*! 
GO TO 70 

911 IF! ICON! ICR ).l E.C) GO TO 901 


COMPUTE CnhERENCE 

CALL RFADMS(9,RI,NSPCT, Jl 1 
CALL REA0MS!9,R1 (NSRCTtl) ,NSPCT, J2) 
CALL READMS!9, SPECT.NPT ,riCHP«-ll 
CON=DELT AT'**2/ ! fc. 283 1 85 30«* WC ON* N8LK. ) 
00 902 r = l,.NSPCT 
XPLOT! I 1 = !I-1) /TMAX 
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OEN3M=,25*(RI ( I )*RI(''lSPrTM ) » 

IF(UtmiM) ^tl2, ^12,41 J 

412 RT ( I )=0. t GOTO 402 

413 RI ( ! I = CABS( SPECK I I/S;JRt(DENOm) )«C0N 

402 COMTIiMUE 

PRINT 403, FP-AMEL, ( I ,XPL0T( I I ,KI ( 1 ) , 1 =1 , NPR INT ( 

PRINT AND PLOT COHE.iLNCE 

403 FORMA r( / /<< ICOME PENCE * 3 A1 0/ / 1) X« I «3X<=F RE 0ULi-,C Y «4X“i-C0HLR ENL E * 3t OX^I 

i 3X*F'-'EgUENC Y»4X«CnHERr.NCE'« )/(!6,2cl3.5,l6,2rl3.3,I6,2ElJ.y,lo,ECl3 
2.3) ) 

E;NC0nE(oC',<>24,PLA4EL ) TP AC K ( K i ) , T PACK ( X 2 ) 

924 FQRMAT(=:'COHEPENCe * A1 0 , 1 X , A 10 , 2 d X I 
NPCPr)7 = NRCR07U 

WRITF(7) PLABEL ,N3PCT , ( XPl ni ( 1 ) , RI ( 1 ) , 1 = 1 ,NSPC T ) 

PRINT 1000, NP.CRD7,PLABEL 
ENCODE! 50,4C4, PLABEL I 

404 FORMAT! <<C0HERt:NCe*41X) 

CALL PLOTNR lYl. ABEL , F R A M E L , 2 3 , XPLOT , R 1 , N S PCT , v , F 1 , F 2 , PL A EL , 1 0 , I E F , 
. 11 ) 

401 IF! ITRA! IC'R ), EQ.C ) GO TO 70 

COMPUTE transfer FUNCTION 

CALL .READMS!9,SPECT,NPT,NCHPf 1) 

IFUTRA! ICP ).LT.0) go TO 403 
CALL REAQMS !9, P, 1 ,NSPCT , J 1 ) 

ENCODE! 6 0,406, PL AH EL ) TRACK !K1 ) ,TKACK(K2 ) 

GO TO 40 7 

405 CALL READMS !9, R I ,NSPCT , J2 ) 

ENCODE ! oO ,408, °LAREL ) T P AC K ! K 1 ) , T PAC K ! .< 2 ) 

4C8 F0RMAT!«TRANSFF.R FUNCTION, TRAYX FOR *A10,3H X ,AiO,dX) 

406 FORMAT (♦TRANSFER FUNCTION, TRAXY FOR ♦A10,3h X ,A1G,6X< 

407 C0N=UELTAT«^2/ ! 6. 2d3 1 85 306^WCDN^NBLK ) 

DO 409 I =1 ,NSPCT 

XRLCT ! I ) =( I-l ) /TPAX 
DPN0M=.5*RI (f) 

IFIOENOM) 414,414,415 

414 RI ! I )=0. -F GO TO 409 

415 RI ( I ) = CA3S ! SPECT ! I )/DENOM ) ♦CON 

409 CONTINUE 

PRINT AND PLOT TRANSFER FUNCTION 

PRINT 410, PLABEL, ( I ,XPLOT( I ) ,R1 m , 1=1 ,NPRIi\iT ) 

410 FOPMAT (//♦1*6A10, / / 5X* I ♦B X TFR EQU t NCY , ♦ 7X ♦TR 3 ! 9X^ I ♦ 3X*F RE 0 

1UENCY^7X4TRA^ )/!I6,2E13.5,l6,2E13.5,l6,2E13.5,16,2E13.5)) 

NRCR07 = NRCR07«-1 

WRITE!?) PLABEL ,NSPCT, !XPLOT( I ) , RK I ) , 1 = l,NSPCT ) 

PRINT 1000, NRCP07, PLABEL 

CALL PL0TN8 (YL ABEL ,F R AM EL , 2 3 , X PL OT , p I, N SPOT, 0, El, F2, PL A BEL, 24, IFF, 
11 ) 

70 CONTINUE 
RETURN 
END 
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